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Abstract 


This report starts with the analysis of one year of Geosat altimeter data starting from the 
orbits computed with the GEM-T2 potential coefficient model and consistent station coordinates 
(Koblinsky et al., 1990). The first stage in the processing followed the general editing procedures 
implemented by Denker and Rapp (1990) when working with GEM-T1 orbits. Additional 
altimeter data, beyond that used by Denker and Rapp, was selected below -63° latitude, in the 
Mediterranean Sea, and in several areas of high frequency signal. The original radial orbit theory 
is due to Engelis. The analysis solved for corrections to the GEM-T2 potential coefficient model, 
coefficients in a degree 10 potential coefficient expansion, and 8 parameters for each of the 76 arcs 
of data analyzed. The data used included the altimeter data, the GEM-T2 potential coefficients with 
its error covariance matrix, and surface gravity data represented by 1° x 1° mean gravity anomalies. 
The root mean square orbit correction was approximately 75 cm with the corrections to the 
potential coefficients corresponding to geoid changes on the order of 1 18 cm. After applying the 
orbit correction terms the adjusted crossover discrepancies were ± 20 cm with a sample point 
residual of ± 1 9 cm. 

The sea surface topography from this solution did not show the slope problem across the 
northern Pacific Ocean that was seen by Denker and Rapp with the GEM-T1 orbits. Variations of 
the sea surface over the one year of data were analyzed by fixing the geopotential model and orbit 
corrections from the one year solution and solving for a monthly sea surface topography 
representation to degree 15. Variations from the annual degree 15 solutions were analyzed in the 
time domain to find signatures at different frequencies, for example annually and seasonally. 
These changes were also studied to detect local variations of the sea surface. 

In a third stage of analysis a combination solution with the GEM-T2 potential coefficients 
and the recent 30' mean gravity anomaly data set was carried out using the same procedure as 
described by Rapp and Pavlis (1990). The global set of adjusted gravity anomalies was used to 
calculate a potential coefficient model to degree 360. The final potential coefficient model was 
formed by taking the coefficients from degree 2 to 50 from the first combination solution with the 
coefficients from degree 51 to 360 of the last solution. The standard deviations of each coefficient 
were computed from the adjustment process and by error propagation. The cumulative geoid 
undulation commission error of the 91 A model to degree 10, 50, and 360 is 5 cm, 25 cm, and 49 
cm, respectively. 

The OSU91 model was tested through orbit predictions and data fitting; through 
comparisons with geoid undulations computed from Doppler and GPS located stations, and with 
comparisons to geoid undulations implied by Geosat altimeter data. In the latter case the root mean 
square difference between the Geosat undulation (after orbit and sea surface topography correction) 
was 34 cm for OSU91 as opposed to ± 53 cm with OSU89B. 
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1. Introduction 


Denker and Rapp (1990) and Denker (1990) have described the analysis of Geosat altimeter 
data for the recovery of an improved geopotential model, sea surface topography and, primarily, 
improved, in the radial direction, Geosat orbits. The research described in these two papers used 
Geosat data with the GEM-T1 based orbits described by Haines et al. (1989; 1990). The GEMT1 
orbits were computed with the GEM-T1 potential coefficient model and the OPNET Doppler 
tracking data. The tracking station positions were those transformed (not adjusted) into the 
reference frame consistent with that used in the development of the GEM-T1 model. The studies 
noted above demonstrated the validity of the orbit improvement process originally developed by 
Engelis (1987a) and modified for the studies noted. However a problem was identified with the 
degree 1 terms of the sea surface topography that were estimated from the solution. Specifically it 
was found that the (1, 1) terms of the spherical harmonic expansion of sea surface topography 
were incompatible with oceanographic information. This specifically was demonstrated by the sea 
surface slope across the northern Pacific Ocean. The coefficients that were found as part of the 
solution implied a slope different from that expected from historical oceanographic information 
(Levitus, 1982). When the SST (1, 1) coefficients were changed to the ones implied by the 
oceanographic data (Engelis, 1987b) the slope problem disappeared. These results suggested that 
the GEM-T1 orbits might have deficiencies that would cause degree (1, 1) problems in the 
determination of SST. The problems could be caused by inaccurate station coordinates (especially 
for the OPNET stations); inaccurate gravity model coefficients; and perhaps unmodeled effects on 
the satellite orbit. 

Haines et al. (1989; 1990) discussed Geosat orbit determination with the GEM-T2 potential 
coefficient model and other improvements over the GEM-T1 orbits.”One improvement was the use 
of Doppler tracking from several TRANET-2 tracking data and the solution of tracking station 
coordinates for these stations. The preliminary results reported by Haines reported the radial orbit 
error with the GEM-T2 analysis to be on the order of 35 cm as opposed to the 85 cm for the GEM- 
T1 analysis. 

The GEM-T2 Geosat orbits were released in early 1990 (Koblinsky et al., 1990) and 
received by us in April 1990. Because of the problems identified with the GEM-T1 orbits we 
decided to test the GEM-T2 orbits with our software to see if improved (specifically with the ( 1 , 1 ) 
coefficients) sea surface topography could be obtained. Several six day arcs were analyzed on a 
preliminary basis where we found the slope problem had been eliminated with the new GEM-T2 
orbits. A decision was then made to process the full year of Geosat data starting from the GEM-T2 
orbits. 


Denker and Rapp (1990) did not include any surface gravity data in their solution. This 
was deliberately done so that results related to the orbit improvement could be emphasized. 
However it was clear that for any modelling effort to be complete, surface gravity, in terms of 
normal equations, should be incorporated into the new solution. The procedures used for doing 
this will be described on a subsequent section. 

Rapp and Pavlis (1990) described a combination solution of the GEM-T2 potential 
coefficient model, surface gravity data, and gravity information derived from Geos-3/Seasat and 
Geosat altimeter data. The model developed there was complete to degree 360 although a rigorous 
adjustment of the data only to degree 50 took place. 

The models (OSU89A and OSU89B) described by Rapp and Pavlis (ibid) used a terrestrial 
gravity data set developed in July 1989 (Kim and Rapp, 1990). Late in 1990 an update of this data 
base was made in which additional gravity data was incorporated. The updated file is described by 
Yi and Rapp (1991). 
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With the above as background it seemed appropriate to combine the new orbits and new 
data together to come up with, first a new gravity model complete to degree 50, a new sea surface 
topography representation, and improved Geosat orbits. This would then be followed by the 
development of a degree 360 potential coefficient model which would be merged with the results of 
the first step. In essence this report discusses an extension of the work described by Denker and 
Rapp (1990), Pavlis (1988) and Rapp and Pavlis (1990) to arrive at improved estimates for 
numerous quantities. 

2. Satellite Altimeter Data Processing 
2.1 Theory 

We start by a brief review of the theoretical models used to relate the altimeter measurement 
to the parameters being sought. We closely follow (Denker and Rapp, 1990) and define a residual 
sea surface height Ah: 

Ah = h c - p - N c - AN 0 - T - co (2.1) 


where: 


h c — computed ellipsoidal height of the satellite based on the a priori ephemeris; 

p = measured and corrected (for environmental factors) distance from satellite to the sea 
surface; 

N c = geoid undulation based on the same geopotential model used for the ephemeris 
generation; 

AN 0 = neglected (or removed) higher frequency geopotential effects; 

% = tidal effects 

co = oceanic effects (waves, etc.) 

The value of Ah will depend on the corrections to the a priori potential coefficients in two parts. 
The first is through the undulation effects (ANg) and the second is through the effect (Ah G ) on the 
ellipsoidal height through the a priori ephemeris. Let Ahi be the ellipsoidal height error caused by 
initial state errors and other effects and let C be the sea surface topography. Then (ibid, eq. (5)) 

Ah = AN G (AQ m , AS /m ) - Ah c (AQ m , AS/J - Ahi + C (2.2) 


where: 


AC/ m , AS* m are the corrections to the fully normalized a priori potential coefficient model 
of degree / and order m; 

The modeling of AN C , and Ah G is described in Engelis (1987a). We continue with a 
spherical harmonic representation of the sea surface topography although problems with the 
representation are discussed in Denker and Rapp (ibid). We write: 
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(2.3) 


C(<t>, A,) = X X (CfcJcosmX + S^ T sinmX) P/ m (sin<f)) 

t=i m=0 

where <j> and X are geocentric latitude and longitude. The value of Ahj was taken the same as used 
by Denker and Rapp (ibid, eq. (7)): 

Ahi = ao + aicostjrt + a 2 sin\j/t + asAtcosyt 
+ a 4 Atsin\j/t + a5Atsin2\[rt + a6At 2 cos\jrt 

+ a7At 2 sin\jtt (2.4) 


V is the frequency associated with the lcy/rev; t is the time from the beginning of the arc being 
processed; and At is the time relative to the middle of the arc. 

Originally the a^ and aj terms were in the Denker/Rapp procedure because such terms can 
represent resonance terms effects which are not modeled in the GEM-T1 model (a 43 order 
resonance was noted by Haines et al., 1989). Since the GEM-T2 model included a 43 order term 
there was a question if the ag and a 7 terms should be retained. Several arcs were analyzed with and 
without the a6 and a 7 terms. The results indicate that the inclusion of these terms gave significant 
improvement (i.e. the root 'mean square residuals were smaller) over the case that excluded the 
coefficients. For all results to be described in this paper the a6 and a-/ terms were retained in the 
Ahj model. 

2.2 Editing and Computational Procedures 

The initial editing procedures for the Geosat analysis were the same as described in Denker 
and Rapp (ibid, p. 13,153) or Denker (ibid, pp. 10-17). The initial editing deleted data over land; 
data where the standard deviation of a linear fit to 10 per second sea height values was greater than 
10 cm; data where the automatic gain control voltage exceeded 37 db; when the attitude was larger 
than 1°3; when the ocean tide correction was larger than ± 1 m; etc. 

After the data was selected in a 6 day arc using the criteria described above a subsequent 
editing was applied with the following criteria: 

1. Data in shallow seas and continental shelves was deleted if the altimeter measurement fell in a 
30' x 30' cell where the depth was smaller than 1000 m. This criteria is the same as used by 
Denker and Rapp (ibid). 

2. An altimeter data point was deleted (initially) if the geoid undulation contribution from degree 
51 to 360 was larger than 3 m when computed from the OSU89B potential coefficient model (Rapp 
and Pavlis; 1990). In the Denker/Rapp solution the 3 m tolerance was applied to the 37 to 360 
contribution so that the new criteria accepts more data. In addition the OSU89B model is more 
accurate than the OSU86F model used by Denker/Rapp so that more appropriate editing results. 

3. An altimeter data point was deleted (initially) if the along track deflection of the vertical 
exceeded 10" as computed from the actual sea surface height data. This criteria is the same as used 
by Denker/Rapp. 

4. Nq data below -64° latitude was deleted because of latitude considerations. Such data was 
deleted by Denker and Rapp (ibid) because of suspected inaccuracies of the OSU86F reference 
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model. Since the OSU89B model is substantially better than 86F, especially in the Antarctic 
regions, we are now able to process data in this region. 

The purposes of criteria 1 was to eliminate data where the tide model might not be as 
accurate as other areas. Criteria 2 and 3 were used to delete data where there was a concern on the 
effect of unmodeled high frequency effects on the coefficients being solved. 

After data was selected with the above criteria sample points were calculated. A sample 
point was evaluated by first fitting the 1 sec residual sea surface height (eq. (2.1)) to a linear fit 
over 20 secs. An iterative outlier rejection criteria with a 3a-limit was used. The sample point 
value was computed from the linear fit at the mid time point of the interval. 

In this computation the N c value was computed from the GEM-T2 potential coefficients 
which are complete to degree 36 with many coefficients extending to degree 50. The initial 
computation of the geoid was for the "zero geoid" (Rapp, 1989) so that an additional correction (to 
be added to the zero geoid) was needed to obtain the "mean geoid" which would be consistent with 
the tidal system used to define the sea surface height. The correction was (ibid, eq. (17)): 

AN = - 0.198 0. sin 2 (J) - jjm (2.5) 

The value of ANo was computed from the OSU89B coefficients from degree 51 to 360 plus 
values for the coefficients not included in GEM-T2. 

Using the editing and selection procedure described in the above paragraphs sample points 
were formed for 76 six day arcs over 22 cycles of the Exact Repeat Mission. The total number of 
sample points was approximately 0.8 million or an average of 10526 points per arc. 
Approximately 20 million one second data were analyzed to arrive at this sample point count. The 
number of normal points will vary from arc to arc depending, primarily, on season or time of year. 
A typical normal point plot from a 17 day Geosat repeat cycle is given in Denker and Rapp (ibid. 
Figure 2). A listing of the Geosat arc number, arc start and stop times and other information as 
reproduced from Koblinsky et. al (1990) is given in Appendix A. 

The above procedures exclude data in several areas that may be of importance to our final 
solutions. Specifically the editing criteria deleted data where the tide was greater than 1 m. Such a 
procedure would delete data where no tide value was present on the GDR. This would 
systematically delete data in some areas where the tide magnitudes are small and the altimeter data 
otherwise useful. The largest such area is the Mediterranean Sea where the tides are small typically 
being on the order of 20 cm although reaching 1 m in some limited areas. In order to bring the 
previously deleted data into the solution, 3 or 4 arcs from each of 4 widely space in time ERM's 
were selected for analysis. The arcs were 4, 5, 6, and 7 (approximately ERM 2); 22, 23, and 24 
(approximately ERM 7); 39, 40, and 41 (approximately ERM 12); and 66, 67, 68, and 69 
(approximately ERM 20). Normal points were computed on the basis of 10 second linear fits in 
contrast to the 20 second interval used previously. This was done to increase the number of 
normal points in this specific region. For these points the ocean tide value was set to zero. For the 
14 arcs a total of 687 normal points were obtained. 

A second inappropriate data exclusion area was the generic area that was deleted because 
the high frequency contribution was significant. Additional consideration recognizing a more 
accurate high degree model (OSU89B) is now available (than used by Denker/Rapp) suggested that 
data be selected in the areas in which such data was previously deleted. Three regions where high 
frequency information was present were chosen as follows: 
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Area One: - 58° £ <}> £ - 14°; 165° <,"K<, 190° 


Area Two: - 60° ^ <|) ^ - 50°; 300° ^ X ^ 340° 

Area Three: 50° £<^60°; 165°<X<210° 

These areas are primarily ones in which trenches can be found. The 14 arcs described earlier 
(ERM 2, 7, 12, and 20) were used to select 701 normal points using the 20 second linear fit 
procedure. 

The standard deviation for each normal point was computed as in Denker and Rapp (ibid, 
p. 13,155): 

SD(Ah) = (0.2 2 + a 2 ) 1/2 m (2.6) 


where 0.2 m is somewhat arbitrary and a is the root mean square misfit (in meters) from the 
normal point computation. 


2.3 The Normal Equation Formation 

The edited data described in the previous section was used to form normal equations that 
will be combined with the GEM-T2 error covariance matrix and the surface gravity normal 
equations (see Section 3). In this formation 76 arcs, over 22 ERMs, were processed. The total 
number of potential coefficient parameters being estimated was 2595. For each arc, 8 parameters 
were estimated leading to 608 arc dependent parameters. 

The parameters associated with the sea surface topography representation were dependent 
on the maximum degree of the expansion as expressed by equation (2.3). Two maximum degrees 
were used in the analysis. In one series of tests the maximum degree was 10 leading to 120 
parameters while in another series of tests the maximum degree was 15 leading to 255 parameters. 
Recall that the zero degree term is deleted from the SST expression effectively forcing the mean 
value to be zero. 

The formation of the normal equations was carried out based on the standard deviations 
computed from eq. (2.6) for all points except for the specially selected points in the Mediterranean 
Sea area and the points from the areas of significant high frequency information. For these points 
the standard deviation computed from eq. (2.6) was multiplied by two. In the case of the 
Mediterranean data this was done to compensate for the use of 10 sec normal points instead of 20 
sec points used in all other areas. In the case of the high frequency content points the 
multiplication by two was used to take into account a higher uncertainty caused by uncertainty in 
the high frequency signal removed from the normal point observation. Other weighting schemes 
could be used or tested but time factors did not allow us to do this. 
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3. Normal Equations for the Geopotential Coefficients Obtained from Surface 
Gravity Data 

In the present study the formation of normal equations for the complete set of geopotential 
coefficients up to harmonic degree 50, from the analysis of terrestrial gravity measurements, 
closely followed the modeling and estimation procedures discussed in detail by Pavlis (1988). 
Therefore, in the following paragraphs only a brief outline of these procedures will be given, and 
the emphasis will be placed on the description of the improved gravity anomaly data which were 
used in this analysis, on certain aspects of the modeling which were re-examined and modified, 
and on the presentation of the results obtained. 

3.1 The OSU October 1990 Gravity Database 

The fundamental terrestrial 1° x 1° mean gravity anomaly dataset used in this study is 
designated "OSU October 1990" (Yi and Rapp, 1991), and represents the latest update of the 
global gravity anomaly database maintained at the Ohio State University. With respect to its 
predecessor (OSU July 1989 — Kim and Rapp, 1990), it is improved by the incorporation of 
improved 944 1° x 1° mean gravity anomalies. Of the 944 newly accepted values, 9 anomalies had 
no previous estimates and 935 anomalies replaced previous values which were primarily 
geophysically predicted. The OSU July 1989 dataset, on the other hand, is substantially better 
than the June 1986 file (Despotakis, 1986) which was originally used for the normal equations 
formed by Pavlis (1988). For example, improved gravity anomaly data for Africa, included in the 
July 1989 dataset, have replaced corresponding values in the June 1986 file, which were identified 
to be contaminated by significant systematic errors (Pavlis, 1988, section 5.3.2). 

The October 1990 database contains in total 50802 1° x 1° mean free-air gravity anomalies. 
Of these, 45932 values originate from actual gravity measurements, while 4870 values are 
estimates obtained from geophysical prediction techniques. In Table 1, statistics related to the 
mean anomalies of the October 1990 dataset are given, while in Figure 1 the geographic 
distribution of the available data is displayed. 

Table 1. Statistics of the 1° x 1° Mean Free-air Anomalies in the OSU October 1990 Database 



Gravity 

Measurements 

Geophysical 

Prediction 

Combined 

dumber of values 

45932 

I* 4870 


Percentage of area 

79.2 

6.6 


Minimum value 

- 270 

-123 

-270 

Maximum value 

303 

127 

303 

Mean value 

-0.5 

-1.0 

-0.5 

RMS value 

27.6 

25.3 

27.4 

RMS standard deviation 

12.0 

17.3 

12.5 


(Gravity anomaly units are mgals; mean and RMS values given above are weighted by the area of 
each 1° x 1° block.) 

3.2 Creation of Mean Anomaly Files Input to the Adjustment 

The analysis made by Pavlis (1988) has demonstrated that the geophysically predicted 
anomalies are in many cases systematically biased with respect to the anomalies that are implied by 
global geopotential models derived from the analysis of satellite perturbations only (ibid, section 
5.3.2). Pavlis and Rapp (1990) have shown that a preferable alternative to the use of geophysically 
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igure 1. Geographic Distribution of the 50802 1° x 1° Mean Free-air Anomalies in the October 1990 Database, “x’’ Identities 
Values Originating From Gravity Measurements (45932) and Geophysically Predicted Values (4870). 



predicted anomalies in global gravity modeling, is the use of anomalies evaluated from the 
combination of lower-degree harmonics obtained from a satellite-alone global gravity model, 
augmented by higher-degree harmonics of topographic/isostatic induced potential. This technique 
was implemented in the development of the OSU89B global geopotential model (Rapp and Pavlis, 
1990), where the coefficients of GEM-T2 (Marsh et al., 1990) up to degree 36 were augmented by 
the topographic/isostatic coefficients of SET3 (Pavlis and Rapp, 1990, p. 373), from degree 37 to 
degree 360, and the resulting "combined" set was used to evaluate 30' x 30' mean anomalies for 
areas occupied by geophysically predicted values, or devoid of any anomaly estimate. 

Based on the results of the above studies, it was decided to adopt the following strategy for 
the preparation of a mean anomaly file to be used as input for the normal equation formation: 


(a) The coefficients of GEM-T2 up to degree 9 are augmented by the coefficients of SET3 from 
degree 10 to degree 50, to form a "combined" set of coefficients Cnm- A global set of 1 ° x 1° mean 
anomalies, Ag" 1 *, are then evaluated by: 


("-’lift 

AOj (rf) n =2 vf / m =- n 


C nm IYU, 


(3.1) 


where the indices (i, j) identify the location of the 1° x 1° block in a two-dimensional array with 
i = 0, 1, ... , 179 and j = 0, 1, ..., 359. (For notation definitions see (Pavlis, 1988)). The rational 
behind the reduction of the "cut-off" degree from 36, used in the development of OSU89B to 9 
used here, is as follows. The normal equations to be produced here are subsequently combined 
with the normal equations that produced the GEM-T2 model. This combination is performed 
under the assumption that GEM-T2 and the geopotential coefficients obtained from the current 
analysis of surface gravimetry, represent two uncorrelated estimates of the true coefficients of the 
field. To account for correlations between these two estimates is (currently) not feasible due to 
computational limitations. It is thus preferable to reduce the higher degree of the harmonics of 
GEM-T2 used in the evaluation of Ag TI , to better comply with the above assumption of zero 
correlation. In addition, the study of Pavlis and Rapp (1990, section 4.2) has indicated that Ag TI 
evaluated using 9 as "cut-off degree are not substantially worse than those evaluated using 36. 

(b) A merging process was performed next whereby a Ag TI value was used to provide the mean 
anomaly estimate for a given 1° x 1° block if: 

— The October 1990 estimate for the block originates from geophysical prediction, or, 

— No estimate is available in the October 1990 database aM the 1° x 1° mean elevation of 
the block is positive. 

In this manner, the 45932 1° x 1° mean anomalies of the October 1990 dataset originating from 
actual measurements (denoted Ag OCT9 °) are maintained in the merged file (denoted SET A), while 
Ag 11 are used to "fill-in" the remaining land areas . SET A contains in total 54048 1° x 1° mean 
anomaly estimates covering 87.3% of the area of the Earth. Of these, 8116 values (8.2% of the 
Earth's area) are Ag TI estimates. The geographic distribution of the data in SET A is shown in 
Figure 2. 

A number of systematic corrections need to be applied to the Ag OCT9 ° anomalies in SET A 
before these can be used in the formation of normal equations. These are (Pavlis, 1988): 


(i) Atmospheric correction 5gA 

(ii) Ellipsoidal corrections £h, £y, £p 

(iii) Second-order vertical gradient of normal gravity correction 8gh2 
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Figure 2. Geographic Distribution of the 54048 1° x 1° Mean Free-air Anomalies in SETA, “x” Identifies Values Originating 
From Gravity Measurements (45932) and “1” Ag^ 1 Values (81 16). 





(iv) Gravity formula transformation 5g r 


The specific formulation used to evaluate 1° x 1° area-mean values of the above correction terms is 
given in detail in (ibid, section 2.3). It needs to be mentioned here that the OSU89B geopotential 
model, complete to degree 180, was used to evaluate 1° x 1° area-mean values of the ellipsoidal 
corrections (denoted IEh, IE-^, and IEp). Also, throughout this analysis the following parameters 
were used to define the geometry and the gravity potential of the reference ellipsoid (ibid, p. 60): 

a = 6378137. m 

1/f = 298.257 

GM = 3986004.36 x 10 8 m 3 /s 2 
co = 7.2921 15 x lO- 5 rad/s 

and the transformation of the October 1990 anomalies, which refer to the GRS 1967 gravity 
formula, to the gravity formula implied by the above constants was performed as explained in 
(ibid, pp. 60-61). Denoting by Ag' the corrected anomaly, one has: 


Ag 


>j 


Agf™ 



(3.2) 


where the total systematic correction 8g s (ibid, equation 4.12): 

(Sgs)ij = [SgA - (IEh + IEy + IEp) + 8gh2 + 8g r ]ij . (3.3) 

The Ag' values in SET A represent surface mean free-air anomalies in the Molodensky 
sense. Their frequency content is not uniform worldwide but depends on factors such as the 
distribution of gravity measurements inside each 1° x 1° block and the averaging process used to 
estimate each mean value. In contrast, the Ag 71 values in SET A are formally interpreted as mean 
free-air anomalies continued to the surface of the reference ellipsoid, and their spectral content 
extends (by definition) only up to harmonic degree 50. Extensive analysis discussed by Pavlis 
(1988, section 5.2.5) has shown that the leakage of power from the higher-frequency component 
of Agjj, to the lower-frequency coefficients being solved for from the incomplete set of discrete 
area-mean values Agjj, can be minimized by removing the higher-frequency content of Agjj (above 
degree 50) prior to the formation of normal equations. The higher-frequency component, 8g 111 ', can 
be evaluated in terms of 1° x 1° area-mean values, using an existing high-degree geopotential model 
such as OSU89B, by: 



Cnn?IY 


li 

rim 


(3.4) 


In the implementation of this procedure the following two aspects need to be considered carefully: 

1. The harmonic coefficients used to evaluate 8g HF must represent as precisely as possible the 
higher-frequency component of the data which will be used in the adjustment. 

2. 8g I1F must be evaluated at the same level at which the mean values Agjj refer (the topographic 
surface of the Earth). 

To comply with the first of the above requirements (since no high-degree expansion was available 
at the time, that included the gravity data from our new source), the following steps were taken: 
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(a) Each of the 944 1° x 1° blocks, containing a value from the new data source, was divided into 
the four corresponding 30' x 30' blocks that it covers. To each of them a 30’ mean anomaly 
estimate was assigned, which is identical with the corresponding 1° x 1° mean value. A total of 
3776 Ztg30' were thus produced. After application of the atmospheric, ellipsoidal, second-order 
gradient and gravity formula corrections and analytical downward continuation to the reference 
ellipsoid (Rapp and Pavlis, 1990), these "split-up" values were merged with the adjusted global 
30' mean anomaly file that produced the OSU89B geopotential model (ibid). In this merging, a 
"split-up" value from China replaced a previous 30' value only if the previous estimate was a "fill- 
in" anomaly or a "split-up" from previously available 1° x 1° estimate. In this manner, 3669 (out of 
the 3776) China 30' mean values were accepted in the merged file. Of these, 3144 values replaced 
"fill-in" data and 525 values replaced previous "split-up" data. 


(b) The resulting global 30' mean anomaly file from the above merging was harmonically analyzed 
as explained by Rapp and Pavlis (ibid, equation 20) to yield^a "modified OSU89B" set of 
coefficients complete to degree and order J60. This set, denoted C^ B , was then used in equation 
(3.4) to evaluate 1° x 1° mean values of 5g\ As it can be seen from (3.4), the high-frequency 
contribution to the data is evaluated on the surface of the reference ellipsoid, while the data Agij, 
refer to the topographic surface of the Earth. To account for this incompatibility, the analytical 
continuation term gi (Wang, 1988) was used in two different ways, which led to two alternative 
files to be used as input for the formation of normal equations: 

Method 1 : Input anomalies referring to the topographic surface of the Earth. 


The global 30' x 30' set of gi values computed by Wang (ibid) on the basis of the TUG87 
mean elevations (Wieser, 1987), was harmonically analyzed (according to a quadrature formula 
similar to equation (20) in Rapp and Pavlis (1990)), to yield _a set of (spherical) harmonic 
coefficients G nm , complete to degree 360. Using the coefficients G nm , two correction terms were 
evaluated, both in terms of 1° x 1° mean values: 


A<Ji (rfr n=2 \lf / m=-n 


GnmlY^n 



1 GM 


360 


AO; (rf) 2 n=51 


2 ( n * 1) 



GnmlYl 


(3.5) 


(3.6) 


The anomalies to be input to the adjustment, denoted Ag (1) , were then defined by: 

A 0 ) h^ 790 + l 5gs ^ ' ( 8g " F ’ s”)ij 

(AgJ 1 - (gi L )ij (3.7) 

depending on their origin (actual measurements or topographic/isostatic values). From the 
definition of the analytical continuation term gi (Wang, 1988), it can be seen that AgO ) refer to the 
topographic surface of the Earth. The anomalies Ag(0 constitute the file designated SET 1 , which 
was one of the two files considered as input to the least-squares adjustment. 
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Method 2 : Input anomalies reduced to the reference ellipsoid. 


Here, the 1° x 1° mean values of the term gj (denoted gi) as computed by Wang (1988) are 
used, and the anomalies to be input to the adjustment are defined by: 

Ag<2>- p8S Cre0 + ( 5 &)ij + <8l)ij-(8s HF )ij 

W? (3.8) 


The anomalies Ag( 2 ) constitute the file designated SET 2, and represent values reduced to the 
surface of the reference ellipsoid. 

In theory, the two alternative treatments of the input anomaly data should yield the same 
result, provided that the observation equations appropriately consider the surface to which the data 
refer. In practice however, the approximations involved in the evaluation of the gi terms (Wang, 
1988), and the errors introduced in the computation of gj L and g* 1 result in small but systematic 
differences as it will be seen later. It should be mentioned here that the harmonic decomposition of 
gi, which is necessary to evaluate gj L and gj H , is not a trivial step since gi has a discontinuity at 
the continental boundary (gi = 0 over the ocean). 

In Table 2 statistics related to the anomalies of SET 1 and SET 2, as well as their 
differences, are given. When forming these files, all Ag^ 1 anomalies were assigned identical 
standard deviation of 20 mgal, based on the accuracy assessment for these values discussed by 
Pavlis and Rapp (1990). Also, the minimum standard deviation for any anomaly regardless of 
source was set to 2 mgal to avoid over optimistic accuracy estimates. The distribution of data in 
both SET 1 and SET 2 is obviously identical to that of SET A given in Figure 2. 

Table 2. Statistics of the 1° x 1° Mean Anomalies Used in the Normal Equations Formed. 



■ESiKHi 



Number of values 

54048 

54048 

54048 

Percentage of area 

87.3 

87.3 

87.3 

Minimum value 

-199.4 

-199.4 

-10.4 

Maximum value 

141.3 

141.2 

9.2 

Mean value 

-0.3 

-0.1 

-0.2 

RMS value 

19.0 

19.2 

0.7 

RMS standard deviation 

13.0 

13.0 

— 


3.3 Estimation of Geopotential Coefficients from Surface Gravity Data 


The anomalies Ag(’) and Ag( 2 ) previously defined lead to the following observation 
equations respectively: 


50 n 

AOi Tjj n =0 Fy/ m=-n 


(3.9) 
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(3.10) 



where v® is the residual associated with the Ag\jp observation (k = 1. 2) and QL represent the 
adjusted geopotential coefficients obtained on the basis of surface gravity data alone (even zonal 
harmonic coefficients are remainders after subtraction of the coefficients of the normal potential). 
The geocentric distance fjj in (3.9) was evaluated as explained in (Pavlis, 1988, section 2.3.4) 
using the TUG87 (Wieser, 1987) 1° x 1° mean elevations to realize the topographic surface of the 
Earth, r? on the other hand is the distance from the geocenter to the point on the ellipsoid at the 
mid latitude of the (i, j) lh block and thus possesses equatorial symmetry. The primes of the 
summations in (3.9) and (3.10) indicate absence of the first-degree terms. The inclusion of the 
zeroeth-degree term is necessitatedjpy the fact that the incomplete set of discrete mean values used 
gives rise to covariances between Cj Q and the rest of the coefficients, which must be taken into 
account (Pavlis, 1988). 

Both observation equations (3.9) and (3. 10) are of the form: 


V = AX - L b 


(3.11) 


Minimization of the weighted norm of the residuals (V T PV) under the condition (3.1 1) yields the 
normal equation system: 

(A t PA)X= A T PL b (3.|2) 

and the least-squares estimate X is: 

X = (ATPA)-l ATpL b . (3.13) 


In the above_A is the design matrix, X is the vector containing C„ m , L b is the vector of 
observations Ag' ; j (k = 1,2) and the weight matrix P is defined by: 


P = 



(3.14) 


with Oq being the a-priori variance of the unit weight (taken to be 1) and ELb 1 * 10 variance- 
covariance matrix of the observations. The a-posteriori variance of unit weight is given by: 

g 2 = Y t PY 

° d.f. (3.15) 

where d.f. are the degrees of freedom, and the variance-covariance matrix of the estimates is: 

= <j?(a t pa)-' 

(3.16) 

For the purpose of combining the normal equations obtained here with corresponding 
normals obtained from the analysis of satellite perturbations, as well as with normals from altimeter 
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measurements, it is critical that s L bP ro P erl y reflects the accuracy of surface gravity data. 
However, the inhomogeneity of the data sources based on which global mean anomaly databases 
are compiled (Kim and Rapp, 1990) makes it difficult to provide realistic estimates of the anomaly 
error variances, let alone error covariances between the mean values. In addition, consideration 
even of simplistic models for the error covariances would make the formation of normal equations 
practically impossible (for degrees of expansion equal to 50 or higher), due to computational 
limitations. Accordingly, following previous experiences (Rapp and Pavlis, 1990), it was decided 
to consider a diagonal Lu> matrix, but modify the original error estimates for the anomalies in an 
attempt to compensate for the neglected error covariances. Denoting by Gjj the standard deviations 
of the anomalies in SET 1 (which are identical to those in SET 2) and by Gq the modified values 
used to form Eu>, the following relationship was imposed: 

max(8, 2 x G°) < of] <min(16, 2 x Gjj) . (3.17) 

This modification yields a ratio 4 : 1 between maximum and minimum weights used in the 
adjustment, and approximately corresponds to the weighting scheme used by Rapp and Pavlis 
(ibid) for 30' x 30' mean anomalies. According to the modification (3.17), the RMS standard 
deviation of the anomalies input to the adjustment is 13.7 mgal, so that the overall accuracy of 
either SET 1 or SET 2 remains practically unchanged (see Table 2). With (and thus P) being 
diagonal, the formation of A T PA and A T PL b can be done efficiently using analytical expressions 
that avoid the use of matrix algebra, as explained by Pavlis (1988, section 4.2.1). 

According to the above, two sets of normal equations were formed using SET 1 and SET 2 
as input data respectively. Both normal equation sets correspond to an expansion complete to 
degree and order 50. From each normal system the corresponding coefficient estimates were 
computed. These are designated VI (from SET 1) and V2 (from SET 2). In both cases 2598 
unknown coefficients are estimated on the basis of 54048 1° x 1° mean anomalies, so that the 
degree of freedom (d.f.) is 51450. In Table 3, statistical information related to the solutions VI 
and V2 is given. 


Table 3. Statistical Information Related to the Gravity Solutions VI, V2 and Cl, C2. 



VI 

V2 

Cl 

C2 

'-l 

Go 

Min Vjj (mgal) 

Max vjj (mgal) 

Mean vjj (mgal) 

RMS vjj (mgal) 

Number of Ivjjl > 7 mgal 

k 

0.302 

-117.8 

187.5 

0.0 

7.3 

11116 

0.303 

-117.7 

187.5 

0.0 

7.3 

11118 

0.354 

-119.5 

195.2 

0.4 

8.2 

14200 

1.078 

0.353 

-119.3 

195.3 

0.2 

8.1 

14207 

1.089 


As it can be seen from Table 3, the solutions VI and V2 are only marginally different (as 
expected). The average percentage difference between them is 4.4%, while the RMS undulation 
and anomaly differences are 1.14 m and 0.73 mgal respectively. In Figure 3 the locations of the 
11118 residuals from V2 which exceed in magnitude 7 mgal are shown. It is clear from this figure 
that V2 fits well the input data over well surveyed (gravimetrically) continental areas (North 
America, Australia, Europe and Africa), while most of the large residuals occur in ocean areas. 
This is primarily due to the incompatibility between the high-frequency component of the surface 
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anomalies over the ocean, with the corresponding component of the altimetry-derived anomalies 
that are used in the evaluation of the "modified OSU89B" coefficients (see also (Pavlis, 1988)). It 
should be emphasized here that the residuals from the solutions employing only surface gravity 
data, represent dominantly a "goodness-of-fit” of the estimated coefficients to the input data. 
Long-wavelength errors that may be present in the surface anomalies cannot be detected without 
the incorporation of superior independent information from satellite-derived normals. 

The statistical information given in Table 3 does not provide any evidence that may be used 
to decide which of the two alternative treatments of the input data (surface values or values reduced 
to the ellipsoid) yields better results. Accordingly, it was decided to compare the results from the 
two methods in greater detail over areas where: 

— High-quality surface gravity data are available 

— Significant variation in elevation is present 

— Independent information, such as undulations obtained from GPS and leveling are available for 
comparisons. 

However, any comparisons between GPS-derived geoidal undulations and undulations computed 
from the solutions VI or V2 may be masked by long-wavelength errors that are present in the 
surface gravity data. For this reason, as well as for the purpose of testing the compatibility of the 
solutions VI and V2 with a satellite-derived model, two preliminary combined solutions were 
performed whereby the coefficients from VI and V2 were combined (in a least-squares sense) with 
the GEM-T2 coefficients (Marsh et aL, 1990). The error variance-covariance matrix 
accompanying each solution was used as weight in order to estimate the "combined" model its a 
weighted average of the two contributing coefficients sets. The resulting "combined" models are 
designated Cl (VI + GEM-T2) and C2 (V2 + GEM - T2). Statistics pertaining to the adjustments 
that produced Cl and C2 are given in Table 3. As it can be seen the RMS residuals from these 
adjustments are only by about 10% higher than the corresponding values obtained when fitting the 
surface gravity data alone. This provides an overall measure of the compatibility between the 
terrestrial and the satellite implied solutions. In Figure 4 , the locations of 14207 residuals from 
C2, exceeding in magnitude 7 mgal are shown. In this figure extended areas in Asia and South 
America are identified, where the terrestrial and satellite-implied anomalies are in disagreement. 
Note that many of these areas cannot be identified in Figure 3, since the residuals from the 
terrestrial-only solution represent primarily a goodness-of-fit to the data as explained before. As 
part of the combination adjustments, calibration factors (see Section 4) were also computed, 
considering the "combined" models versus GEM-T2 as a subset solution. These are given in Table 
3 and their values indicate that the weighting scheme used for the surface anomalies, yields 
satisfactory results. The average percentage difference between the solutions Cl and C2 is 2.7%, 
while the RMS undulation and anomaly differences are 0.08 m and 0.36 mgal respectively. 

Using the harmonic coefficients of the solutions Cl and C2, geoidal undulations were 
computed next according to: 



where the notation definitions are given in (Rapp and Pavlis, 1990, pp. 21899-21900). The 
differences AN 12 = N(C1) - N(C2) over the areas of Europe and North America are shown in 
Figures 5 and 6 respectively. From these figures it can be seen that AN 12 are highly correlated 
with elevation. The signatures of the Alps in Europe and of the Rocky Mountains and Sierra 
Nevada in North America are clearly identifiable. 
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Figure 6. Undulation Differences N(C1) - N(C2) in North America. Point Values Computed on 1° x 1° Grid and 
Plotted with Contour Interval Cl = 10 cm. 





In the area of Europe 12 GPS stations forming part of the European GPS traverse (Torge et 
al„ 1989) were selected in the area of interest, in order to compare the geoidal undulations derived 
from GPS positioning and leveling to the corresponding values obtained from the solutions Cl and 
C2. For this purpose the contribution to the gravimetric undulation from degree 51 to 360 was 
computed using the coefficients of the model described in Section 6. The results from these 
comparisons are summarized in Table 4. 

Table 4. Undulation Comparisons at Selected GPS Stations in Europe. ANj = Ngps - N(C1), 

AN 2 = N gps - N(C2). 


Station 

Number 

<pO 

X(°) 

mmmssm 

Bill 

(m) 

Ngps 

(m) 
AN | 

(m) 

an 2 

1 

50.9227 

9.7604 

250.190 

47.64 

0.16 

0.23 

2 

50.5089 

9.6814 

336.324 

48.25 

0.28 

0.31 

3 

50.2988 

10.4542 

336.478 

47.55 

0.28 

0.28 

4 

50.0894 

10.9823 

318.680 

47.04 

-0.01 

-0.03 

5 

49.7782 

11.0287 

299.920 

47.09 

0.19 

0.14 

6 

49.3477 

11.0215 

392.799 

47.07 

0.33 

0.23 

7 

49.0361 

11.1316 

604.813 

47.21 

0.72 

0.58 

8 

48.6780 

11.5877 

439.395 

46.56 

0.63 

0.45 

9 

48.3248 

11.5709 

540.154 

46.23 

0.69 

0.47 

10 

48.0417 

11,6324 

613.217 

46.22 

0.85 

0.60 

11 

47.7803 

11.7235 

804.996 

46.94 

0.91 

0.63 

12 

47.4901 

11.2521 

954.553 

49.13 

1.33 

1.00 




Mean Difference 

0.53 

0.4 1 




Standard Deviation Diff. 

0.37 

0.26 


From Table 4 it is evident that N(C2) is in better overall agreement with Ngps tlian N(C1). 
It is recognized here that the sample of the 12 GPS stations used for the comparisons is too small 
to support a definitive argument as to wether the surface values or the values continued to the 
ellipsoid provide a better modeling of the terrestrial anomaly data. The undulation differences 
between the two alternative solutions, being on the order of 20 to 40 cm over mountainous areas, 
require very accurate independently derived undulations so that meaningful comparisons can be 
made. Nevertheless, based on the limited evidence presented here, it was decided that the normal 
equations formed considering the anomalies continued to the ellipsoid (solution V2) are to be 
preferred. These normal equations were subsequently used in combination with GEM-T2 and the 
normals obtained from satellite altimetry, to provide the solutions discussed in section 4. 

The combination of two normal equation sets obtained from different data is straight 
forward in case the two sets were formed using the same approximate values for the unknown 
parameters (Pavlis, 1988, p. 68). If this is not the case, one of the two sets needs to be 
"translated" to the approximate values of the second. To illustrate the principle let: 


NX = U (3-19) 

be the normals obtained here from surface gravity, where: 

N = A T PA ; U = A T PL b (3.20) 

These refer to the ellipsoidal even zonal coefficients as approximate values, so that (3.19) can be 
written as: 
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(3.21) 


n(Xtot - Xell) = u 

where Xtot are the adjusted gravitational harmonic coefficients and Xell the vector containing the 
values of the even zonal coefficients of the ellipsoidal gravitational potential (and zeroes for the rest 
of the coefficients present in Xtot)- If X x2 contains the values of the coefficients in GEM-T2, 
equation (3.21) can be written as: 

n[(Xtot - Xto) + (Xt2 - Xell)] = U (3.22) 

or: 

n(Xtot - X-n) = U - N(Xt2 - Xell) (3.23) 

so that, to refer the normals (3.21) to the GEM-T2 approximate values the vector U needs to be 
"translated" by - N(Xt2 - Xell)- The principle (obviously) applies to any change of approximate 
values. 
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4.0 The Initial Combination Solutions to Obtain a Potential Coefficient Model to 
Degree 50 

The general least squares adjustment procedure, with a priori parameter weights, was used 
in the estimation of the unknowns. The equations are described in Rapp and Pavlis (1990, p. 21 , 
889). For this discussion we define the normal equation form as follows: 

N = A t PA (4- 1 ) 

where A is the design matrix and P is the weight matrix assumed here to be a diagonal matrix. We 
have two types of normal equations: one from the altimeter data, N A ; the other from surface 
gravity data, Nq. The solution vector (actually the correction vector to the a priori parameter 
values) is: 

V X = - (N a + N g + Pt 2 + Px) _1 (U A + Ug) ( 4 - 2 ) 

where Pn is the inverse of the error covariance of the GEM-T2 potential coefficient model, Px arc 
the a priori weights on the selected parameters of the adjustment. U A and Ug are the misclosure 
vectors. 

In our adjustment the a priori weights were used in two cases: 1) sea surface topography 
coefficients and 2) potential coefficients not included in GEM-T2. In the sea surface topography 
coefficients we need to fix (see discussion in Denker and Rapp, 1990, p. 13,158) the degree 1, 0 
term. This term can not be separated from the 1-cy/rev orbit correction term and therefore we fix 
the value at an oceanographic estimate. This values is 0.1297 m based on the harmonic analysis of 
the Levitus sea surface topography by Engelis (1987b, Table 1, p. 28, Ocean Solution to Degree 
10). The initial value of the other coefficients were taken as zero with standard deviations (to 
degree 10) based on the root mean square coefficient implied by the ocean solution. For SST 
solutions to degree 15 the standard deviation for coefficients from degree 1 1 to 15 was taken to be 
the same as that at degree 10. The standard deviations used are given in Table 5. Also included in 
this table are the standard deviations implied by the SST signal model given by Nerem et al (1990a, 
eq. 20). With the exception of degree 1, the standard deviations used in this paper are roughly 
50% larger than used by Nerem et al. This simply implies that we put less of a constraint on our 
SST estimates than was described in the Nerem et al study. 

Table 5. A Priori Standard Deviations for Each SST Harmonic Coefficient (Unit = cm) 


Degree 

This Paper 

Nerem 

1 

12.8 

19.0 

2 

22.4* 

7.88 

3 

10.5 

4.71 

4 

6.41 

3.27 

5 

6.22 

2.46 

6 

6.36 

1.95 

7 

4.35 

1.61 

8 

2.66 

1.35 

9 

2.29 

1.17 

10 

1.31 

1.02 

11 

1.31 


12 

1.31 


13 

1.31 


14 

1.31 


15 

1.31 



* based on a degree variance of 2500 cm 2 . 
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A priori estimates for the C 21 and S 21 potential coefficients were also included. These 
coefficients were not incorporated in the altimeter observation equations but they were in the 
gravity anomaly observation equations. Consequently the gravity normal equations have such 
coefficients and correlations do exist between these coefficients and others in the solution. To 
force the final C 2 .i and S^i coefficients to be zero a high weight was assigned to force the zero a 
priori value. 

This analysis used one year of Geosat data represented by 22 ERMs. Since each ERM 
represents the same geographic region, the use of 22 ERMs would cause a disproportionate fit of 
the model to the altimeter data. Consequently a down weighting of the altimeter data is needed to 
assume a balanced solution. The down weighting problem was discussed by Denker and Rapp 
(1990, p. 13,156). 

Initial combination solutions were made with down weighting factors of 1/24 and 1/96 
using the initial VI surface gravity normal equation set as described in Section 3. Solutions were 
made with sea surface topography expansions to degree 10 and 15. Using the improved V2 
gravity normal equations sea surface topography solutions to degree 10 and 15 were also made 
with the 1/24 and 1/96 down weighting of the altimeter normal equations. Table 6 defines the 
various solutions made for this study. 

Table 6. Designation of One Year Solutions 


Solution 

SST(max) 

Alt. Weights 

Surface Normals 

FYS10.W24GV 

10 

1/24 

VI 

FYS10.W96GV 

10 

1/96 

VI 

FYS15.W24GV 

15 

1/24 

VI 

FYS10.W24GW 

10 

1/24 

V2 

FYS10.W96GW 

10 

1/96 

V2 

FYS15.W96GW 

15 

1/96 

V2 


As noted in Section 3 the solutions corresponding to the VI normal equation set were 
found slightly poorer in several tests than the V2 solution. Therefore, in subsequent analysis our 
results will primarily refer to the V2 solution. 

We first compare the effect of the altimeter down weighting factor on the corrections to 
residual sea surface heights, geoid undulation corrections, crossover discrepancies, etc. Results 
for six arcs out of the 76 used are shown in Table 7. All values are given for the solution of sea 
surface topography to degree 10. In this table C is the root mean square sea surface height implied 
by the harmonic coefficients and Ah (residual) is the root mean square residual of the altimeter 
measurement observation equation, eq. 2.2. 

Table 7 can be compared to Table 4 in Denker and Rapp (1990). The comparison will 
show that the value of Ah G + Ah[ has been reduced in this new solution by about 30% and AN G has 
been reduced by 51%. The sea surface topography, rms residuals, and rms crossover 
discrepancies are all comparable in both solutions. We also see from Table 7 that the results are 
insensitive to the 1/24 or 1/96 weighting scheme. We do see a small (1 cm) increase in the Ah 
(residual) and crossover discrepancy when going from the 1/24 to 1/96 weighting. The increase is 
expected; the small rms is welcome. 
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Table 7. Root Mean Square Value of Corrections to the Various Terms in the Error Model, and 
Crossover Discrepancies Based on the Adjusted Orbits and the One Year Solution. 
Units are cm. 


Arc 

Solution 

AhQ + Ah i 

an g 

C 

Ah Model 

AhResidual 

AdjCrossover Diff. 



71 

119 

64 

143 

18 

19 

H 

KSrtKvKRSsvB 

71 

118 

63 

142 

19 

20 


Ri lilHKE HEM 

82 

115 


149 

19 

21 


Ikl It 

81 

114 

BS1B 

149 

19 

22 

3 

S10.W24GW 


Ql 

66 

152 

18 

21 


S10.W96GW 

Bail 1 

mam 

65 

151 

19 

23 

22 


70 

119 


rnmnam 

19 

21 


BdMfrMMSl 

70 

118 



20 

22 

23 


73 

117 

82 

Hi i 


21 



72 

116 

81 

MUM 

Wmmm 

22 

24 

wxbwmmm 

77 

117 

83 

169 

18 

20 



76 

116 

82 

169 

19 

21 


The next comparison was to compare the heights at the one day overlaps in the six day 
arcs. The results are shown in Table 8 for the degree 10 SST case. 

Table 8. Mean and RMS Value of Satellite Height Differences at One Day Overlapping Arc 
Segments. Units are cm. 


| Difference 


Solution 

Mean 

RMS 


S10.W24GW 

2 

5 


S10.W96GW 

2 

5 

2/3 

S10.W24GW 

-4 

5 


S10.W96GW 

-4 

5 

22/23 

S10.W24GW 

-3 

4 


S10.W96GW 

-3 

4 

hh 

S10.W24GW 

4 

6 

BBMBaiaHMM 

S10.W96GW 

4 

6 


This table gives results that indicate there is no sensitivity to the altimeter weighting scheme 
in the orbit overlap comparisons. 


Next we calculated the calibration factors for several solutions with respect to the GEM-T2 
model and the standard deviations of the coefficients implied by the error covariance matrix. The 
calibration factor was introduced by Lerch et al. (1988) and used by various investigators (e.g. 
Rapp and Pavlis, 1990, p. 21, 898) to calibrate various models. Let AF n be the root mean square 
(rms) coefficient difference, at degree n, between two potential coefficient models. Specifically: 


AF n = 


ACI 


AS nm 


life 


Lm=0 


N r 


(4.3) 


where N n is the number of coefficients at degree n. The rms coefficient error at degree n, for the 
solution case would be: 
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(4.4) 


CJn = 


£ °fc ^) +o L) 


-m=0 


Nr, 


1 /2 


with a similar expression, for the GEM-T2 model. We define a term £ as follows: 
£l = o£ 2 - 

For a given degree the calibration factor is: 


(4.5) 


k n — AFn/^n 


(4.6) 


The average calibration factor is: 


k = 



50 

XN n k n 


n=2 


(4.7) 


Values of k are given in Table 9. 


Table 9. Calibration Factors (k) 


Solution 

Value 

FYS10.W24GW 

1.16 

FYS15.W24GW 

1.15 

FYS10.W96GW 

1.14 


The values in Table 9 indicate the low sensitivity of the calibration factor to the altimeter weighting 
factor. Assuming the weighting procedures used in the solution are appropriate, the value of k = 
1.14 indicates the standard deviations implied by the GEM-T2 error covariance matrix may be 
slightly optimistic. (Recall that the ideal value of k would be one.) 

In order to further investigate the impact of the different altimeter weighting procedures the 
standard errors of the geoid undulations and sea surface topography were computed. This was 
done by calculating the propagated error in the parameters of the solution to the undulation error 
and the sea surface topography error on a 3°6 grid. The root mean square error was calculated for 
land and ocean areas. An ocean area was defined to be where -70° < <f> < 70° and h < -200 m, 
using a 30' x 30' elevation file. The results are given in Table 10 for geoid undulation and Table 
16 in Section 5 for sea surface topography. Also shown are global averages computed from the 
gridded values and the error degree variances. These two values will be similar but not identical 
since the use of the error degree variance approach neglects the error correlation of the coefficients. 

Table 10. Standard Error of Geoid Undulation Based on Selected Degree 50 Combination 
Solutions. Units are cm. 


Solution 

Ocean 

Land 

Global 

Through Degree Variances 

FYS10.W24GW 

9.6 

36.5 

23.2 

21.9 

FYS10.W96GW 

14.2 

39.4 

26.2 

24.8 


From Table 10 we see that the undulation error in the ocean areas is 26% of the undulation 
error in the land areas. The use of the 1/96 weighting procedure increases the undulation error 
from 9.6 cm (1/24 case) to 14.2 cm. A smaller percentage change is seen for the undulation error 
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on land. Also note that the error computed through the degree variances is slightly (6%) smaller 
than when the more rigorous approach is used. 

A comparison was made of the geoid undulations implied by the 1/24 and 1/96 wci gluing 
solutions. The differences are shown in Table 11. 

Table 11. Comparison of Geoid Undulations Implied by the FYS10.W96GW and 
FYS10.W24GW Solution to Degree 50. Units are in cm. 



Ocean 

Land 

Global 

RMS Dif. 

6.4 

22.3 

14 

Minimum Dif. 

-61.4 

-92.6 

-92.6 

Maximum Dif. 

43.0 

94.7 

94.7 


The RMS differences seen in this table are quite consistent with the standard errors of the solutions 
seen from Table 10. Of interest is how the weight changes on the altimeter data effect the 
undulations on land. 

At this point several external evaluations were carried out to help in the selection of the 
preferred altimeter weighting scheme. A complete discussion of external comparisons will be 
made in Section 8 but two pertinent comparisons will be addressed here. In these comparisons, 
the degree 50 potential coefficient model (W96GW or W24GW) was augmented by the coefficients 
of the OSU89B model from degree 51 to 360. 

The first comparison was with geoid undulations computed at the GPS stations described 
in Rapp and Pavlis (ibid, p. 21,900-21,901). Table 12 shows the standard deviation of the 
undulation difference (GPS derived minus model) for the two different weighting procedures. 

Table 12. Standard Deviation of GPS/Leveling Implied Geoid Undulation Minus Model Geoid 
Undulation for Two Altimeter Weighting Procedures. Units are cm. 


[ Mode! 


Traverse 

FYS10.W24GW 

FYS10.W96GW l 

Europe 

37 

34 

Canada 

37 

37 

Australia 

39 

35 

Scandinavia 

35 

32 

Tennessee 

24 

21 


The values given in Table 12 imply the 1/96 weighting gives slightly better comparisons at the 
GPS sites. Examination of geoid undulation difference maps showed differences that could reach 
95 cm although the rms undulation difference between the two solutions was 14 cm. 

The next test involved the comparison of the geoid undulation implied by a Geosat ERM 
with the undulation implied by the geopotential model. The comparison is generally described by 
Rapp and Pavlis (ibid, Section 4.4), In the comparison to be described here the satellite orbits and 
sea surface topography to degree 10 were implied by a preliminary solution with a 1/24 weight. 
Table 13 shows these comparisons for 768286 points on approximately ERM 7. 
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Table 13. Comparison of Potential Model and Altimeter Implied Geoid Undulations 


| Model 



FYS10.W24GW 

FYS10.W96GW 

Std. Dev. 

33.2 cm 

34.0 cm 

No. > 1.5 m 

4149 

4386 

No. > 2.0 m 

1326 

1437 


The numerical results indicate a slightly poorer fit when the 1/96 weighting was used. This is 
clearly expected but the only very slight deterioration indicates the 1/96 can be used without any 
significant deterioration of the altimeter fits. 

Taking into account the information given in the past several tables the 1/96 weighting 
procedure was selected to be that used for the OSU91A solution. The main reason for the selection 
was related to the increased undulation discrepancies at the GPS stations when using the 1/24 
procedure. This, coupled with the only slight altimeter fit deterioration when the 1/96 weights 
were used, led to the 1/96 weighting procedure for the final solution. 

Figure 7 shows the geoid undulation commission error implied by the FYS10.W96GW 
solution. This map can be compared to Figure 8 in Denker and Rapp (ibid) or Figure 1 1 of Rapp 
and Pavlis (ibid). However it should be noted that the Figure 7 map reflects the errors in the 
complete degree 50 model while the other figures were based on coefficients to degree 36 that were 
in the GEM-T1 model. We see from this figure the significantly improved accuracy in the ocean 
areas as compared to the land. This, of course, was also apparent, from Table 10. The geoid 
accuracy in the well surveyed (gravimetrically) areas (U.S., Europe, Australia) is on the order of 
30 cm while the accuracy is 50 to 60 cm in the poorer areas. Of obvious note is the smaller errors 
in the polar regions despite the fact that data coverage is sparse. A reasonable explanation for why 
this error is lower is needed. The accuracy by degree, and cumulatively, will be discussed in 
Section 6. We can repeat, from Table 10, that the global undulation accuracy of the model (1/96) 
is 26 cm which may be contrasted with 33 cm for the OSU89B model to degree 50 (Rapp and 
Pavlis, ibid, Table 11). 

The bias (or a Q ) term in equation (2.4) was examined for all 76 arcs. For the (1/24) 
weighting the bias was 62.0 ±2.1 cm while for the (1/96) weighting it was 62.3 ± 2.1 cm 
indicating no essential bias difference in the weighting schemes. Adopting the (1/96) bias value 
would imply an equatorial radius of 6378137.00 m minus 0.623 m or 63781 36.38 m. The bias 
values were individually examined and plotted in Figure 8 starting with the first arc. The bias 
value is plotted at a time associated with the middle of the arc. A cyclic variation with an amplitude 
of about 2 cm is seen. The data was spectrally analyzed with the spectrum shown in Figure 9. 
There are two clear lines; one at 2 cycle/year and the other at 22 cycles/year. The first line is 
associated with a semi-annual variation while the second line is associated with the length (17.05 
days) of a Geosat ERM. The key point is to understand why there should be a semi-annual 
variation in the bias term. Is the variation connected with seasonal pressure changes on the oceans; 
is it associated with changes in the tropospheric corrections; do seasonal winds play a role; are 
there tide errors with a semi-annual signal; etc. Additional study is needed to explain the bias 
variation. 

This Section has examined the first stage in our development of a geopotential model to 
degree 50, a sea surface topography model to degree 10, and orbit correction parameters for 76 
Geosat arcs covering the first year of the ERM. Different weighting procedures were tested, 
evaluated, and selected for further use. The root mean square rms residual for the altimeter 
observation was typically ± 19 cm with the rms crossover discrepancy being on the order of ± 22 
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Figure 8. Geosat Bias Term (ao) for 76 Arcs Based on the FYS10.W96GW Solution. 



Figure 9. Spectrum of Geosat Bias Terms for the 76 Arcs Used in the FYS 10.W96GW Solution. 
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cm. This is a considerable improvement from the 49 cm crossover discrepancies associated with 
the original GEM-T2 orbits (Haines et al, 1990, p. 2884), In the next section we consider details 
of the sea surface topography solutions. 



5.0 The Sea Surface Topography Representations 

5.1 The One Year Average 

As noted in the previous section a number of solutions for sea surface topography were 
made. Expansions to degree 10 and degree 15 were made with our final altimeter weighting 
scheme of 1/96. For each solution the degree (1, 0) term was fixed at the value implied by the 
Engelis harmonic analysis of the Levitus SST information. The harmonic coefficients and their 
standard deviation are given in Table 14 for the degree 10 and Table 15 for the degree 15 solution. 

In Figures 10 and 1 1 the sea surface topography for the degree 10 and degree 15 solution 
are given. Figure 10 can be compared to Figure 6 Denker and Rapp (1990). The essential 
difference between the two is that the solution of this report shows a reasonable slope across the 
northern Pacific Ocean, while the Denker/Rapp initial solution did not. In order to achieve a 
reasonable slope, Denker/Rapp were forced to fix the degree (1, 1) harmonics at oceanographically 
implied values. In the current solution no such fixing is required. This does imply that the GEM- 
T2 orbits, with a consistent gravity model and geocentric station set, yields better long wavelength 
oceanographic information than the GEM-T1 orbits. 

Another diffeience between the 10, 10 solution of this paper and the Denker/Rapp solution 
is the high near 0° latitude, and 210° longitude. This high was not apparent in Denker/Rapp but it 
has been seen in the solutions (PGS3853) of Nerem et al. (1990). 

The geostrophic currents implied by the degree 10 solution are shown in Figure 12. This 
map can be compared to Plate 1 of Denker and Rapp (ibid) where similarities exist but numerous 
differences can be seen. For example, the California current is more reasonably represented on the 
newer map than the earlier version. There is a complex circulation pattern around latitude 0° and 
longitude 210°. The solution is finding it difficult to separate the North Equatorial Current from the 
Equatorial Counter Current. This is not unreasonable since a degree 10 solution corresponds to a 
linear resolution of 1998 km. Other major currents (e.q. Kuroshio Current, Mozambique Current, 
Gulf Stream (in a broad sense). South Equatorial Current, Brazil Current, Antarctic Circumpolar 
Current, etc.) are evident from the flow map. However there is no sign of the Agulhas Current off 
the south tip of Africa. The degree 15 flow vector map shows more detail than seen in Figure 1 2. 
However these details are fragmented and not quite believable. For the record it is shown as 
Figure 13. 

We have examined in greater detail sea surface topography implied by the degree 10 
solution in the Mediterranean Sea. Figure 14 shows the SST with a 10 cm contour interval. This 
figure shows a decrease of the SST in a northeast direction. The maximum change in SST in the 
Mediterranean Sea, from Africa (<t> = 31°, X = 18°) to the north end of the Agean Sea is 98 cm. 
From the entrance to the Mediterranean to the north end of the Agean Sea, the SST difference is 94 
cm. 


Figure 14 may be compared to Figure 31 of Lisitzin (1974) which displays the mean sea 
level in the Mediterranean Sea. The slope pattern follows quite well that shown in Figure 14. For 
example, the mean sea level change from the entrance to the Mediterranean to the Agean Sea is 80 
cm (Agean is lower). This compares well with the 94 cm found in this study. 

We have evaluated sea surface topography estimates from several other harmonic models in 
the Mediterranean Sea. None show patterns resembling the Lisitzin result. This is because there 
SST solutions did not incorporate any data form the Mediterranean region in their solution. 


31 


Table 14. Sea Surface Topography Harmonic Coefficients and Their Standard Deviation Implied 
by the First Year of die ERM, Degree 10 Model. Units are meters. 


n 

m 

C 

s 

CT ( C ) 

a ( s ) 

l 

0 

0 * 1297 


0.0010 


2 

0 

- 0.4933 


0.0341 


3 

0 

0.1667 


0.0292 


4 

0 

- 0.0582 


0.0277 


5 

0 

0.1045 


0.0257 


6 

0 

0.1818 


0.0229 


7 

0 

- 0.0312 


0.0189 


8 

0 

0.0466 


0.0150 


9 

0 

- 0.0374 


0.0119 


10 

0 

0.0101 


0.0090 


1 

1 

- 0.1550 

- 0.0978 

0.0268 

0.0333 

2 

1 

- 0.0628 

0.0366 

0.0173 

0.0409 

3 

1 

- 0.0415 

- 0.0308 

0.0179 

0.0392 

4 

1 

0.0231 

0.0015 

0.0185 

0.0305 

5 

1 

- 0.0056 

0.0149 

0.0160 

0.0273 

6 

1 

0.0294 

0.0119 

0.0164 

0.0241 

7 

1 

0.0912 

0.0290 

0.0142 

0.0205 

8 

1 

- 0.0559 

- 0.0389 

0.0119 

0.0160 

9 

1 

0.0217 

0.0429 

0.0117 

0.0140 

10 

1 

0.0372 

- 0.0224 

0.0084 

0.0095 

2 

2 

0.0106 

0.0324 

0.0200 

0.0181 

3 

2 

0.0329 

- 0.0511 

0.0281 

0.0209 

4 

2 

- 0.0104 

- 0.0400 

0.0293 

0.0216 

5 

2 

0.0107 

- 0.0024 

0.0255 

0.0206 

6 

2 

- 0.0008 

0.0188 

0.0208 

0.0181 

7 

2 

0.0285 

0.0302 

0.0172 

0.0158 

8 

2 

- 0.0028 

0.0145 

0.0144 

0.0134 

9 

2 

0.0000 

- 0.0267 

0.0123 

0.0118 

10 

2 

- 0.0207 

0.0017 

0.0095 

0.0093 

3 

3 

- 0.0174 

- 0.0269 

0.0165 

0.0177 

4 

3 

0.0282 

- 0.0107 

0.0191 

0.0203 

5 

3 

0.0077 

- 0.0266 

0.0204 

0.0221 

6 

3 

0.0156 

- 0.0718 

0.0189 

0.0207 

7 

3 

0.0255 

- 0.0356 

0.0153 

0.0163 

8 

3 

0.0201 

0,0064 

0.0127 

0.0132 

9 

3 

- 0.0343 

0.0282 

0.0116 

0.0117 

10 

3 

0.0074 

0.0150 

0.0084 

0.0082 

4 

4 

- 0.0048 

- 0.0265 

0.0134 

0.0157 

5 

4 

- 0.0337 

0.0180 

0.0170 

0.0162 

6 

4 

- 0.0014 

- 0.0094 

0.0166 

0.0170 

7 

4 

0.0179 

- 0.0229 

0.0151 

0.0159 

8 

4 

0.0319 

- 0,0162 

0.0124 

0.0127 

9 

4 

0.0079 

- 0.0032 

0.0106 

0.0106 

10 

4 

0.0108 

- 0.0168 

0.0078 

0.0077 

5 

5 

0.0039 

0.0050 

0.0129 

0.0135 

6 

5 

0.0173 

0.0029 

0.0139 

0.0141 

7 

5 

0.0282 

0.0024 

0.0135 

0.0135 

8 

5 

0.0012 

0.0054 

0.0117 

0.0120 

9 

5 

0.0221 

- 0.0172 

0.0101 

0.0102 

10 

5 

- 0.0117 

0.0070 

0.0074 

0.0075 

6 

6 

0.0109 

- 0.0179 

0.0124 

0.0100 

7 

6 

0.0244 

0.0161 

0.0111 

0.0114 

8 

6 

0.0081 

- 0.0050 

0.0101 

0.0103 

9 

6 

- 0.0049 

- 0.0056 

0.0092 

0.0090 

10 

6 

- 0.0092 

0.0121 

0.0068 

0.0066 

7 

7 

0.0205 

0.0056 

0.0099 

0.0086 

8 

7 

- 0.0060 

0.0053 

0.0097 

0.0090 

9 

7 

0.0018 

0.0098 

0.0083 

0.0081 

10 

7 

- 0.0013 

0.0065 

0.0065 

0.0062 
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8 

8 

- 0.0030 

0.0151 

0*0085 

0.0081 

9 

8 

0.0024 

0.0075 

0.0079 

0.0082 

10 

8 

- 0.0180 

0.0103 

0*0060 

0.0060 

9 

9 

0.0098 

0.0108 

0.0074 

0.0074 

10 

9 

- 0.0094 

0.0165 

0.0063 

0.0060 

10 

10 

0.0135 

- 0.0012 

0.0064 

0.0069 


The standard error of the sea surface topography is shown for the degree 10 solution in 
Figure 15. This map is comparable to Figure 9 of Denker/Rapp. From this figure we see that the 
standard deviation is fairly uniform over the oceans, increasing as one approaches land areas. The 
standard error of the SST has been computed in the ocean areas (-70° < ^ < 70°); h < -200 m) for 
three solutions developed for this report. There values are shown in Table 16. 

Table 16. Standard Error of Sea Surface Topography in the Ocean Areas. Units are cm. 


Solution 

Standard Deviation 

Through Deg Variances 

FYS10.W24GW 

±7.0 

14.9 

FYS10.W96GW 

8.4 

16.6 

FYS15.W96GW 

11.6 

20.6 


From Table 16 we see that the standard deviation of the sea surface topography in the ocean 
areas increases from 7.0 cm with the 1/24 weighting to 8.4 cm with the 1/96 weighting (for the 
degree 10 case). These values are approximately 50% of that found through the error degree 
varianced. Recall that the error degree variances reflect the higher standard deviations in the land 
areas. The error for the 15, 15 solution is higher because twice as many coefficients are included 
in the degree 15 model. 

It is usual to plot the power spectrum of the sea surface topography and compare it to 
corresponding error (SST and geoid) spectra. Denker and Rapp (1990) have pointed out serious 
limitations of this process since the SST is defined only over the oceans. Nevertheless, for 
continuity purposes we show in Figure 16, the square root of the power spectra of the two sea 
surface topography solutions, as well as the errors, by degree, of these models and the geoid 
undulation errors that are part of solution leading to the degree 10 SST model. 

From Figure 16 we see that up to degree 6, the spectrum from the degree 15 solution is less 
than that of the degree 10 solution. For example, at degree 10, the signal from the degree 10 
solution is 6.9 cm while it is (at the same degree) 5.0 cm from the degree 15 solution. This 
suggests some spectral leakage from the higher harmonics into the lower harmonics. After degree 
10, the signal from the degree 15 solution is nearly flat. 

The degree errors of the degree 15 solution are always greater than the degree 10 solution. 
At degree 10, it is 30% higher. The geoid undulation error is always smaller than the signal 
suggesting that SST determination out to degree 15 could be made. However the discussion in 
Denker and Rapp (1990, p. 13,157) notes the high negative correlation between the SST 
coefficients and die potential coefficients (geoid information) between degree 10 and 15. This case 
would still exist for the solutions described in this report. 

We next examine the second degree zonal harmonic of sea surface topography. In principle 
this harmonic is independent of the tidal (non-tidal, zero, mean) reference system, if the sea surface 
and geoid are consistently treated. Some papers, however, have published C 2 ,o values in non- 
consistent systems. However, knowing the system, one can convert to a common system. Some 
of the conversion procedures are described in Rapp (1989). Other solutions are described in 
Nerem et al (1990, p. 3167) and Marsh et al. (1990, p. 13,143). The SST coefficients given in 
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Table 15. Sea Surface Topography Harmonic Coefficients and Their Standard Deviation Implied 
by the First Year of the ERM. Degree 15 Model. Units are meters. 


T1 

m 

C 

s 

o ( c ) 

a ( s ) 

1 

0 

0.1297 


0.0010 


2 

0 

- 0,5173 


0.0375 


3 

0 

0.1899 


0,0328 


4 

0 

- 0.0761 


0.0300 


5 

0 

0.1070 


0.0288 


6 

0 

0.1659 


0.0252 


7 

0 

- 0.0188 


0.0222 


8 

0 

0.0186 


0.0175 


9 

0 

- 0.0351 


0.0154 


10 

0 

0,0099 


0.0108 


11 

0 

0.0029 


0.0101 


12 

0 

- 0.0125 


0.0090 


13 

0 

- 0.0009 


0.0084 


14 

0 

- 0,0008 


0.0082 


15 

0 

0.0007 


0.0077 


1 

1 

- 0.1540 

- 0.0913 

0.0286 

0,0350 

2 

1 

- 0.0585 

0.0207 

0.0190 

0.0428 

3 

1 

- 0.0455 

- 0.0410 

0.0199 

0.0411 

4 

1 

0.0320 

0.0148 

0.0203 

0.0319 

5 

1 

- 0.0108 

0.0055 

0.0178 

0.0288 

6 

1 

0.0243 

0.0265 

0.0179 

0.0258 

7 

1 

0.0685 

0.0251 

0.0176 

0.0232 

8 

1 

- 0.0380 

- 0.0356 

0,0146 

0.0180 

9 

1 

0.0206 

0.0151 

0.0143 

0.0166 

10 

1 

0.0259 

- 0.0020 

0.0103 

0.0111 

11 

1 

- 0.0108 

- 0.0022 

0.0101 

0.0108 

12 

1 

0.0121 

0.0044 

0.0092 

0.0099 

13 

1 

- 0.0034 

- 0.0186 

0.0093 

0.0096 

14 

1 

- 0.0060 

0.0218 

0.0084 

0,0090 

15 

1 

0.0113 

0.0000 

0.0079 

0.0084 

2 

2 

0.0160 

0.0296 

0.0223 

0,0213 

3 

2 

0.0335 

- 0,0540 

0.0296 

0.0227 

4 

2 

0.0036 

- 0.0381 

0.0306 

0.0232 

5 

2 

0.0060 

0.0009 

0.0267 

0.0222 

6 

2 

- 0.0066 

0.0160 

0.0221 

0.0196 

7 

2 

0.0335 

0.0184 

0.0193 

0.0180 

8 

2 

0.0110 

0.0012 

0,0163 

0.0156 

9 

2 

- 0.0207 

- 0.0170 

0.0149 

0.0147 

10 

2 

- 0.0016 

- 0.0006 

0.0106 

0.0105 

11 

2 

0.0182 

- 0.0072 

0.0104 

0.0104 

12 

2 

0.0002 

- 0.0079 

0.0092 

0.0093 

13 

2 

- 0.0105 

0.0033 

0.0093 

0.0093 

14 

2 

0,0218 

- 0.0080 

0.0090 

0.0091 

15 

2 

- 0.0152 

0.0086 

0.0091 

0.0091 

3 

3 

- 0.0143 

- 0.0185 

0.0190 

0.0203 

4 

3 

0,0169 

- 0.0248 

0.0209 

0.0217 

5 

3 

0.0129 

- 0.0247 

0.0220 

0.0234 

6 

3 

0.0096 

- 0.0575 

0.0204 

0,0221 

7 

3 

0.0088 

- 0.0419 

0.0173 

0.0180 

8 

3 

0.0206 

- 0.0050 

0.0148 

0.0150 

9 

3 

- 0.0221 

0.0240 

0.0136 

0.0134 

10 

3 

- 0.0027 

0.0084 

0.0104 

0.0103 

11 

3 

- 0.0025 

- 0.0084 

0,0099 

0.0098 

12 

3 

0.0020 

0.0021 

0.0094 

0.0092 

13 

3 

0.0022 

0.0056 

0.0094 

0.0093 

14 

3 

- 0.0016 

- 0.0070 

0,0089 

0.0087 

15 

3 

0.0119 

- 0.0046 

0.0086 

0.0085 

4 

4 

- 0.0118 

- 0.0250 

0.0153 

0.0186 

5 

4 

- 0.0270 

0.0105 

0.0186 

0.0179 
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6 4 

7 4 

8 4 

9 4 

10 4 

11 4 

12 4 

13 4 

14 4 

15 4 

5 5 

6 5 

7 5 

8 5 

9 5 

10 5 

11 5 

12 5 

13 5 

14 5 

15 5 

6 6 

7 6 

8 6 

9 6 

10 6 

11 6 

12 6 

13 6 

14 6 

15 6 

7 7 

8 7 

9 7 

10 7 

11 7 

12 7 

13 7 

14 7 

15 7 

8 8 

9 8 

10 8 

11 8 

12 8 

13 8 

14 8 

15 8 

9 9 

10 9 

11 9 

12 9 

13 9 

14 9 

15 9 

10 10 

11 10 

12 10 

13 10 

14 10 


0.0045 

0.0165 

0.0272 

- 0.0012 

0.0091 

0.0035 

0.0041 

- 0.0006 

0.0018 

0.0006 

0.0076 

0.0163 

0.0323 

0.0065 

0.0216 

- 0.0118 

- 0.0021 

- 0.0080 

0.0015 

0.0042 

0.0075 

0.0054 

0.0303 

0.0032 

- 0.0056 

- 0.0073 

- 0.0059 

0.0142 

- 0.0018 

0.0082 

0.0064 

0.0128 

0.0106 

0.0011 

- 0.0076 

0.0011 

- 0.0046 

- 0.0012 

- 0.0150 

- 0.0120 

- 0.0028 

0.0096 

- 0.0190 

- 0.0019 

0.0063 

0.0016 

- 0.0093 

0.0017 

0.0116 

- 0.0045 

0.0015 

- 0.0058 

0.0072 

- 0.0083 

0.0022 

0.0126 

- 0.0029 

0.0104 

- 0.0052 

0.0033 


- 0.0078 

- 0.0090 

- 0.0126 

0.0030 

- 0.0117 

0.0074 

0.0006 

0.0099 

- 0.0107 

- 0.0003 

0.0049 

0.0062 

0.0013 

0.0167 

- 0.0171 

- 0.0097 

- 0.0123 

0.0091 

- 0.0017 

0.0035 

0.0008 

- 0.0117 

0.0054 

- 0.0027 

0.0016 

0.0034 

- 0.0063 

- 0.0036 

- 0.0034 

- 0.0123 

0.0117 

- 0.0024 

0.0062 

0.0004 

0.0009 

0.0010 

0.0036 

- 0.0066 

0.0049 

- 0.0098 

0.0139 

0.0111 

0.0179 

0.0031 

0.0054 

0.0034 

- 0.0027 

0.0110 

0.0081 

0.0181 

0.0031 

0.0123 

- 0.0061 

0.0012 

- 0.0121 

- 0.0069 

0.0069 

- 0.0039 

0.0063 

0.0131 


0.0178 

0.0167 

0.0142 

0.0126 

0.0099 

0.0096 

0.0092 

0.0088 

0.0086 

0.0087 

0.0155 

0.0160 

0.0154 

0.0134 

0.0122 

0.0095 

0.0092 

0.0088 

0.0090 

0.0087 

0.0081 

0.0154 

0.0136 

0.0123 

0.0113 

0.0091 

0.0088 

0.0084 

0.0084 

0.0079 

0.0080 

0.0135 

0.0123 

0.0108 

0.0087 

0.0085 

0.0079 

0.0081 

0.0080 

0.0076 

0.0120 

0.0102 

0.0084 

0.0078 

0.0076 

0.0076 

0.0073 

0.0071 

0.0104 

0.0086 

0.0080 

0.0071 

0.0071 

0.0067 

0.0069 

0.0084 

0.0074 

0.0069 

0.0064 

0.0062 


0.0189 

0.0175 

0.0144 

0.0126 

0.0099 

0.0095 

0.0093 

0.0088 

0.0086 

0.0087 

0.0159 

0.0160 

0.0152 

0.0137 

0.0122 

0.0095 

0.0093 

0.0087 

0.0090 

0.0088 

0.0081 

0.0127 

0.0136 

0.0121 

0.0107 

0.0090 

0.0087 

0.0083 

0.0084 

0.0079 

0.0080 

0.0120 

0.0112 

0.0106 

0.0084 

0.0084 

0.0079 

0.0080 

0.0079 

0.0075 

0.0110 

0.0103 

0.0081 

0.0076 

0.0074 

0.0075 

0.0072 

0.0070 

0.0104 

0.0082 

0.0079 

0.0070 

0.0071 

0.0066 

0.0069 

0.0089 

0.0078 

0.0070 

0.0063 

0.0062 
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15 

10 

0.0015 

11 

11 

0.0019 

12 

11 

- 0.0072 

13 

11 

0.0043 

14 

11 

0.0029 

15 

11 

- 0.0071 

12 

12 

0.0047 

13 

12 

- 0.0013 

14 

12 

0.0010 

15 

12 

- 0.0032 

13 

13 

0.0074 

14 

13 

0.0015 

15 

13 

0.0059 

14 

14 

0.0034 

15 

14 

0.0006 

15 

15 

0.0140 


0.0086 

0.0060 

0.0060 

- 0.0092 

0.0081 

0.0081 

0.0105 

0.0070 

0.0072 

- 0.0085 

0.0069 

0.0069 

0.0027 

0.0059 

0.0059 

- 0.0072 

0.0057 

0.0057 

- 0.0035 

0.0071 

0.0071 

0.0095 

0.0064 

0.0063 

- 0.0013 

0.0059 

0.0059 

0.0100 

0.0051 

0.0051 

- 0.0022 

0.0064 

0,0066 

0.0020 

0.0060 

0,0062 

- 0.0061 

0.0051 

0.0052 

0.0001 

0.0061 

0.0062 

- 0.0048 

0.0054 

0.0053 

0.0004 

0.0060 

0.0061 


Table 14 have been derived from a consistent tidal system and no correction term is needed. Table 
17 summarizes some determinations of C 2 ,o for SST. 

Table 17. The Second Degree Fully Normalized Harmonic Coefficient of Sea Surface 
Topography. Units are cm. 


Solution 

Data 

Value 

Mather (1978) 

Geos-3 

-43 

Mather(1978) 

Oceanographic 

-46 

Engelis (1987) 

Oceanographic 

-28 

Ohio State (1988) 

Seasat (17 days) 

-47 

Marsh et al (1990) 

Seasat 

-44 

Nerem et al (1990) 

Geosat (3/17 days) 

-33 

Denker/Rapp( 1 990) 

Geosat (one year) 

-39 

Putney et al (1991) 

Geosat (3 months) 

-34 

TEG-2 (1991) 

Geosat 

-40 

This paper 

Geosat (one year) 

-49 


From Table 17 we see there is a wide range of values for this coefficient. The value obtained in 
this paper is consistent with the ocean results of Mather (1978) and the Ohio State 1988 Seasat 
(17day) solution carried out by Knudsen (1988, private communication). The value is smaller for 
the Geosat solution that has used less data than the new Ohio State solution. 


5.2 Temporal Variations of the Ocean Surface 

The sea surface heights developed from the GEM-T2 improved orbits enables us to study 
the time variations in the ocean surface. For this paper we will examine the variations in two ways 
that are described in the next two sections. The study of large scale time variations has been 
described by Denker (1990, Section 4.4) and Nerem et al. (1990b). 

5.2.1 Temporal Variations Through Data Related to the Harmonic Analysis of Sea 

Surface Topography 

The sea surface topography defined so far represents the average over the first year of the 
ERM mission, November 1986 through October 1987. It is possible to study time variations in 
this year by creating SST expansions for specific time periods using the altimeter data from the 
time period. This expansion can be referred to the one year surface by removing the coefficients 
defined by the one year mean. Then sea surface height differences can be calculated. This 
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LONGITUDE 



Figure 11. Sea Surface Topography Implied by Degree 15 Solution From One Year (1986/87) of Geosat Data. 
Contour Interval is 10 cm. 
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Figure 12. Geostrophic Flow Patterns Implied by the Degree 10 Sea Surface Topography Model. 
Data Calculated On A 5° Grid. 
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Figure 13. Geostrophic Flow Patterns Implied by the Degree 15 Sea Surface Topography Model 
Data Calculated On A 5° Grid 
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Figure 





SQUARE ROOT OF DEGREE VARIANCE (CM) 


100 



HARMONIC DEGREE 


Figure 16. Square Root of Signal and Error Degree Variances for Degree 10 
and Degree 15 SST Models, and Geoid Error. 
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procedure was used by Denker (1990, Section 4.4) who examined differences with 2 month 
averages. He found the differences to be small, approximately 4 cm over the oceans. 

In this study a modification of the procedure developed by Denker was made. First, the 
orbit correction parameters and the adjusted gravity field to degree 50 was held fixed using the (10, 
10) SST model. Then a SST model to degree 15 was computed on a monthly basis starting from 
November 6, 1986. The degree 15 field was chosen to see if time changes at the higher resolution 
(vs degree 10) could be seen. It was recognized that the coefficients between degree 1 and 15 
could be contaminated by constant (in time) geoid error but the time variations would be 
meaningful. In doing this the mean sea surface topography was defined by averaging 12 monthly 
solutions to degree 15. The residual fields were then generated by removing the mean field from 
the monthly field in the spectral domain. The actual values of the residual sea surface topography 
were then calculated on a monthly basis. In Figure 17 the residual SST is shown for the month 
(starting on the 7 th ) of April 1987. The most significant feature is the -20 cm depression near the 


Average Sea Level 

140 - 



Figure 18. Geosat-derived Sea Level, Averaged in Zonal Bands Between 130° and 165°E. Each 

Band Contains Approximately 80 Independent Time Series Expressed as Anomalies 
Relative to the Annual Mean, April 1985-1986. Curves are Offset by 20 cm Relative 
to Each Other. (From Cheney and Miller (1990, p. 2982)). 
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Bands. Units are cm. 


46 


equator at longitude 240°. This residual is weaker (-10 cm) in February March, and stronger April. 
Starting in June 1987 there is no significant deviation from the annual mean in this area. This 
deviation is probably related to the El Nino current in this time period as will be discussed shortly. 
Variations in other geographic regions are on the order of 5 cm, as found by Denker (ibid). Since 
the expansion of SST is only to degree 15 (corresponding to a resolution of 1332 km) mesoscale 
variations can not be seen with this procedure. However large scale variations can be quantified as 
will be discussed in the next section. 

We next examine the sea level changes implied by the degree 15 analysis to see if one can 
detect the 1986-1987 El Nino event. In this study we first considered the results of Cheney and 
Miller (1990) and Miller and Cheney (1990). Figure 6 in Cheney and Miller, shown in this report 
as Figure 18, shows sea surface height residuals with respect to the annual mean of April 1985- 
April 1986. This figure is based on "approximately 80 independent, time series. . .". The data was 
averaged in zonal bands between 130° and 165°E. 

In order to see if our results showed similar signals the relative degree 15 SST was 
evaluated on a 1° x 1° grid. The data was then averaged in the same bands as used by Cheney and 
Miller. In our case this meant an average of 160 points except for the lowest latitude band in which 
234 points were averaged. The results shown in Figure 19 are plotted on a scale similar to that 
shown in Figure 18. Comparison of Figures 18 and 19 indicates a good, but not perfect 
agreement. The positive residual at the two northern latitude bands during May-July 1987 is quite 
consistent in the two figures. The negative anomaly during June-August 1987 for the southern 
most latitude band agrees well in both figures. But our results show a positive anomaly in 
November 1986, a feature not seen in Cheney and Miller (ibid). The change in the residual from 
November 1986 to August 1987 is 20 cm from Figure 18 and is 14 cm from Figure 19. 

One should be cautious in interpreting the positive and negative residuals because the 
absolute zero in the Cheney and Miller analysis is arbitrary. The point here is that temporal 
changes of the ocean surface can be detected through a degree 15 SST representation. The changes 
show a significant agreement with other analysis carried out in a significantly different way. 

Another test carried out was suggested by the results plotted in Fig. 4a of Cheney and 
Miller (1990). This plot shows the variation of sea level in a 8° x 1° cell in the Pacific near Ponape. 
To compute our results the 12 sea surface topography degree 15 representations were evaluated in 
the 8° x 1° cell whose northwest corner was 7° in latitude and 154° in longitude. The 1° x 1° grid 
values were averaged. The average was then removed from the 12 values. The residual plot is 
shown in Figure 20. 



1986 1987 

Figure 20. Sea level Variation Over One Year Based on A Degree 15 Representation of SST. 
Cell's northwest comer is 7° (<t>), 154° ( X ). 
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The diagram shows a low in December (86) with a change of 15 cm to a peak in February (87). 
The Cheney/Miller result shows the low in January with the peak in March. The change is 20 cm. 
Our solution shows a near linear decrease of 13 cm from February to August while the 
Cheney/Miller change is 20 cm from March (87) to September (87). The times appear to be offset 
by a month and the amplitude of the changes somewhat less in our analyses. However the general 
agreement is quite encouraging. 

5.2.2 Temporal Variations Through Spectral Analysis 

As described in previous sections a solution was obtained in which the sea surlace 
topography, the geopotential and the orbit error were simultaneously determined by processing the 
one year Geosat altimeter data. Based on this solution the time variation of the sea surface 
topography will be studied. 

The geoid is assumed to have no time variation in some time period. The assumption may 
not be exactly correct, but it is believed that it is a very good approximation for studying the time 
variation of the sea surface topography and the current geoid determination accuracy. Let sea 
surface topography £ be represented in a series involving time: 

C(<p, X, t) = Yj fem(t)cosm\ + s/ m (t) • sinmA.] P /m (sincp) 

l,m (5.1) 

The coefficients of the spherical harmonics c* m and s/ m are functions of the time. The non time 
dependent version of (5.1) is eq. (2.3). SST may change continuously in time. The frequencies 
of the SST time variation range from zero to infinity. In practice, the SST can only be sampled at a 
time interval for a limited time period. Therefore only the frequencies of a limited band can be 
detected. 

By fixing the geopotential and orbits correction from the one year solution, the SST can be 
expanded into a series for every month. Approximately 6 six days arcs of Geosat data were used 
to form a monthly SST solution: The 12 monthly solutions can be written as 

C(<p, X, t k ) = X fem(tk)cosrnX + s/ m (tk)sinmX.] P/ m (sin<p) 

t,m (5.2) 

k = 0, 1, ... ,11 

where k is the number of the month starting from November 8, 1986. 

Basically, the information of the sea level change during the time period November 1986 to 
October 1987 is included in the 12 monthly SST solutions. Assume that we can fit the SST to the 
following model: 

£(< ?,X,x) = t; o ((p,X) + iiq,X)i (5.3) 

where £ denotes the sea level change rate. The £ and £ are in the form: 


Co(<P> = X ( c^ m cosrnX + s /m sinrnX) P/ m (sincp) 




(5.4) 
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c(<p> aJ = X (c/mcosm^. + s* m sintnX) P/ m (sinq>) 

4m 

and the coefficients c^ m , c/ m and sj m , S/ m are determined by fitting a straight line to the 1 2 sets of 
coefficients c <rn (t k ) and s* m (t k ) in the least squares sense. That is: 

_ o 

c /iri(fk) = C /m + C/ m t k + Rck 

k = 0, , 1 1 

S/m(tk) = sj m + S, m t k + Rsk (5.5) 

O O J_ 

where Rck, Rsk are the residuals of the misfit. The c^ m , c/ m and s^n, s* m will be determined by 
li 

^ R^ k = Minimum 
k=0 


li 

V Rj k = Minimum 

k=o (5.6) 


If the standard errors of the coefficients c/ m (t k ) and s/ m (t k ) are available, condition (5.6) 
becomes 


S 


k=0 


o2 

— 2k- = Minimum 

®c/mk 


n r 2 

y — = Minimum 
k=0 Os/mk 


(5.7) 


The condition (5.6) is a special case of (5.7). It was assumed that the standard error of the 
coefficients are equal to one. The solution with condition (5.7) can be found in many books, e.g., 
Press et al., (1986, p. 504). 


In order to study the energy distribution of the time variation of the sea surface topography 
corresponding to frequencies, eq. (5.2) is written in a discrete Fourier-series: 


C(<p, X, k) = 



(A^mCosrnA, + B n/m sinrnX) P* m (sin(p) e 2 * ,kn /N 

k = 0, .... 1 1 


(5.8) 


where i = V- I. N is the number of the monthly SST solutions and it is assumed to be even. The 
coefficients An/ m and Bn/m are time independent and defined by: 
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An/m-X^‘k) e ' J ' ik " /N 

k=o 

fo (5,9) 

The discrete Fourier transformation and its inverse can be found in many books, e.g., Brigham 
(1988, Chapter 6). The frequency fk is defined by: 

f k = _l_ k = 0, 11 (5.10) 

N * A 


where A is the time interval and for our case that is 1 month. 


Based on the sampling theorem, the Nyquist frequency is given by 


f N = 


1 

2 • A 


(5.11) 


which is 1 eycle/2 month. Therefore the frequency band ranges from 0 to 1 cycle/2 month for our 
SST time variation studies. The physical meaning of the frequencies is: 

k = 0 annual average of SST 
k « 1 annual SST variation 
k * 2 semi-annual SST variation 
k «= 3 4 month S ST variation 
k =? 4 seasonal SST variation 
k m 5 2,4 month SST variation 
k = 6 2 month SST variation 


Eq. (5.8) can also be written in a cosine series. By setting 


Annn “ aijjf m + ibn^ m 

Bfi/jn = + ibn/jfl » (5.1 2) 

we then have: 

(A n r m cosmX + Bn/mSinmX.) 

*= (anrmCosrnX + a n * m sinmX)+ i(bn/ m cosmX + tw m sinmX) ( 5 . 1 3 ) 

Inserting (5,13) into (5.8), we get 

c(<p, X, k) * 1 X X kta + ifcj 

■ n=0 /,m (5.14) 

with 
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0Wm = ( a n^mCostnX + a^mSinmX) P/ m (sin<p) 
Pn/m = (bn/ m cosmX + tw m sinmX) P/m(sin(p) 


(5.15) 


If we write e 2 ™ 1 ^ = cos27tnk/N + isin27tnK/N, and note that the £ is a real function, so that eq. 
(5. 14) becomes: 


C(<P> *>. k ) = jj- X X (“n^mcos 2nkn. + p^sinlska) 

^ n=0 i y m 


(5.16) 


The Fourier coefficients of a real function satisfies the following relationship (ibid, p. 397): 


F N -n = (F n )* (5-17) 

where F n is the Fourier coefficient of a real function, and the asterisk denotes the complex 
conjugate. If we write 


Fn — ttn/m + iPn^m > (5.18) 

then we have from eq. (5. 17) 

a n/m = a (N-n)/m 

Pn/m = *P(N-n)/m (5.19) 

In addition we have 

2rckn 2rck(N - n) 
cos N — cos N 

. 2rckn . 2ttk(N - n) 
sin— jq— = -sm ^ 

By using eqs. (5.19) and (5.20), (5.16) becomes 


C(<P> k ) = ^I a otm +|ZX 5 n' (on/mCOS^ant . p^sin^^j 

/,m n=l /,m 

with 

g 10.5 if n = N/2 
n jl if n= 1, .... N/2 - 1 

If we set 

®n/m = ^n/mCOSOln 

Pn/m = ^n/mSinCOn 


(5.20) 

( 5 . 21 ) 

(5.22) 

(5.23) 
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with 


Qtfm * ®n4n Pn/m 
Pn/m 


cOn = tan' 


ttn/m 


(5.24) 


eq. (5.21) then can be written as a cosine series: 

C(<P> K k) = jjf 2, « 0 /m + £ iL S 5 n Q n /mCOs(2mli + C0 n ) 
W 4m 1 n=l 4m 


(5.25) 


The frequencies contained in (5.25) are positive and range from 0 to the Nyquist frequency fN. 
Eq. (5.25) gives the SST at the time t = k month, k = 0, ... 11. If SST needs to be computed at 
the arbitrary time between Nov. 1986 and Oct. 1987, eq. (5.25) can be easily written as (c.f. 
Cizek, 1986, eq. (5.36b)): 


C(<P> X,, t) = X a o4n + ^ X &n^n/m * 1 + ’ 

N 4 m ^ n=l 4m 


(5.26) 


and SST can be evaluated at the time t by summing the above series. If N is a large number, a 
more efficient computation procedure using FFT is described in (ibid, p. 94). 

The power spectrum of the SST is defined as 

f p((p, X, f n ) = A} £ X 5nQn/mCOSpSjp. + «„) , 

N k=l 4m 

= n = 1, ... , N/2 

N 2 S, 


p{(p, X, f 0 ) — ^o^m n — 0 

N z 4m 


(5.27) 


Based on Parseval's theorem, the sum of the power spectrum of the SST is the square mean of the 
SST: 


S p(*. x. O-iS 

n=0 W k=0 

From eq. (5.28) the mean square of the SST time variation is given by: 

W k=0 
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(5.29) 



where £ denote the mean value of the SST and it is given by 


_ N-i 

c(<p> ^)4I c(<pA,t k ) 

n k=0 

/,m 


(5.30) 


The estimation of the power spectrum of the SST is based on the means of the periodogram 
(Brigham, 1988 p. 367). The statistical meaning of the power spectrum is discussed in (ibid) and 
Elliott et al. (1982). A brief and comprehensive discussion can also be found in Press et al. (1986, 
p. 422). These discussions are interesting and important for interpretation the power spectrum 
defined by eqs. (5.27) and (5.29) correctly but it exceeds the scope of this report, and we simply 
interprete the quantities defined by eq. (5.27) as the power of the SST corresponding to different 
frequencies. 

It is worthwhile to mention that Nerem et al. (1990b) modeled the coefficients of the sea 
surface topography by: 


C/m(0 “ C/m + C/mt (A n /jnCOS(Dit + B n ^ m sinC0jt) 

i=l 

(5.31) 

s/m(t) = s^n + s /m t +X ( 4 , m cosc 0 it + b^ m sincDit) 
i=l 

(5.32) 


o o 

where c#m and s^m are the coefficients of the mean sea surface topography topography. The 
c/ m ands/ m are the coefficients of the sea level change rate. The series are the periodic time 
variation of the SST. 


If the coefficients are determined by using the least-squares technique, the model can fit the 
actual SST well in the meaning of best fit. But the interpretation of the amplitude or power 
spectrum must be carefully made. 

Two years of Geosat ERMs data were used for the SST time variation studies (ibid). The 
frequencies in the ERMs data range from zero to 1 cycle/34.10 days. Only selected frequencies 
(eqs. (5.31), (5.32)) were included in the model. The frequencies which were not modeled in 
(5.31) and (5.32) would be aliased into the model and mixed with the selected frequencies. The 
amplitude of the selected frequencies would be falsified by unmodeled frequencies. Because of the 
aliasing problem it is not meaningful to talk about the annual or semi-annual SST variation and so 
on. 


Another consideration is the separation of the coefficients of the secular SST change 
c/ m ands/ m from the coefficients of the Fourier-series in (5.31) and (5.32). The secular SST 
changes c/ m t and s/ m t are not orthogonal to the base function of the Fourier-series, the cosine and 
sine. It will be of interest to look at the correlation between them in the least-squares solutions. 
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In this section the secular SST changes were solved by using eq. (5.5). The periodic SST 
variation was contained in the terms R c k and R s k which were considered as misfit. One can argue 
that we have the similar situation as described in the last paragraph. The difference is that only the 
secular SST changes were solved by using eq. (5.5). We feel a little safer for the reason that the 
periodic SST variation may effect the straight line fit very little. 

In the following we will give the results of the secular SST change during the period Nov. 
1986 to Oct. 1987 and the temporal variations of the SST through spectral analysis. Figure 21 
shows the yearly sea surface topography change rate computed from eq. (5.7). The SST changed 
considerably in some regions. For example, in the west Pacific ocean a - 1 5 cm/yr rate ol the SS 1 
is found between the longitude 140° to 180° around the equator. This decrease may be due to the 
El Nino phenomenon (Cheney and Miller, 1990). The SST rate was about -15 cm in the Atlantic 
Ocean between latitude 0° to 20°. A small SST change is found in the region of the Brazil and 
Falkland currents. Globally the RMS values of the SST change during the period Nov. 1986 to 
Oct. 1987 is 8.6 cm. The global mean, based on a 1° x 1° grid, of the SST change (this is 
equivalent to the global sea level change) was found to be -0.2 mm in the ocean which is defined as 

f- 70° < (p < 70° 
ocean = { 

[h£-200m , 

where h is depth of the ocean. The very small mean value means that almost no change in the 
oceanwide mean sea level has been detected by our solutions. 

The spectral analyses of the SST time variation enable us to study the SST variation in the 
frequency domain. The variation is decomposed into the components of the annual, semi-annual 
and seasonal variation and so on (cf. eq. (5.25)). The power spectrum of the components of the 
SST variation is half of the square of the amplitude of the Fourier-series in (5.25). It is more 
convenient in the geometric sence to use the amplitude rather than the power spectrum defined by 
eq. (5.27). Table 18 gives the RMS value of the amplitude of the SST time variation 
corresponding to different time periods: 

Table 18. RMS Values of the Amplitude Corresponding to Frequency fk- 


k 

RMS Value 

1 

5.8 cm 

2 

4.2 

3 

2.5 

4 

2.0 

5 

1.7 

6 

0.8 


Table 18 shows clearly that the amplitude becomes smaller, when the frequency gets higher. No 
strong seasonal variation (k = 4) was found. 

Figure 22 and 23 show the amplitude of the annual (k = 1) and seasonal (k = 4) SST 
variation. A relative strong (9 cm) annual variation appears in the west Pacific ocean and around 
the equator. Another relatively strong variation occurs in the Atlantic ocean above the equator. 
The seasonal variation is small. In the western Pacific 2 to 4 cm SST variation can be found. 
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The time variation of the SST was studied on a global basis. We consider the SST as a 
function of the geographic location (<|>, X) and the time (t) and sampled it at the monthly time 
interval. The sampling interval can be smaller, but it is limited by the repeat time of Geosat which 
is 17.05 days. The SST was expanded up to degree and order 15 so that the resolution (~ 1300 
km) is somewhat too broad to identify the detailed SST variation, e.g., the variations due to 
narrow currents. Such variations can be examined by working directly with the altimeter data and 
with the improved orbits. The use of the harmonic analysis will be limited, for now, to broad 
feature changes. These results are preliminary and could be expanded if additional (in time) Geosat 
data were analyzed. 
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Figure 23. Amplitude of the Seasonal Variation (k = 4) of SST Derived from One Year of Geosat Data Based on a Spherical 
Harmonic Expansion To Degree 15. Contour interval = 2 cm. 




6.0 The Estimation of the Potential Coefficients Above Degree 50. 

The previous discussion has described how a potential coefficient model to degree 50 has 
been selected on the basis of a combination of the GEM-T2 potential coefficients, surface gravity 
normal equations, and one year of Geosat altimeter data. We now incorporate these coefficients in 
a development leading to a complete set of potential coefficients to degree 360. The general 
procedure followed is based on the analysis described by Rapp and Pavlis (1990). 

The first step in this process was the updating of the global 30' x 30' mean anomaly file 
used in the OSU89 model development. The first updating was to incorporate the new data in the 
30' data file. The first step was to split the new 1° x 1° mean anomalies into 4 identical 30' x 30' 
values. A total of 3776 30' values were estimated from the original 944 values of anomaly Source 
97. The following corrections were then applied to the 30' values: ellipsoidal corrections; gi 
corrections; transformation to the GEM-T2 implied gravity formula; atmospheric corrections, and 
second order gravity formula corrections. The details of these corrections are given in Section 2.2 
of Rapp and Pavlis (ibid, p. 21,887). 

The corrected 30' file was then merged with the 30' input file used for the development of . 
the OSU89B degree 360 mode! The merging criteria were that a new value would replace a 
previous (OSU89B) value only if the previous value was a “fill-in” value, or it was the result of the 
split up of a previously available 1° x 1° value. More specifically the previous value had to belong 
to either the SET 3 or SET 5 case of Table 2 in Rapp and Pavlis (ibid, p. 21 ,896). The standard 
deviations assigned to the new 30' defined anomalies was equal to that of the value being replaced. 
The code assigned to these anomalies was 4097. The newly created file contained 3669 out of the 
3776 (by split up) 30' anomalies. Of the 3669 values, 3144 replaced fill-ins while 525 values 
replaced previous 30' values from 1° values (see section 2). 

The next step in the updating process was done to reflect a change in philosophy in the 
calculation of fill-in gravity anomalies from GEM-T2 and topographic/isostatic information. In the 
OSU89B development the fill-in anomalies (i.e. anomalies in areas where no data was available) 
were computed from the GEM-T2 model to degree 36 with contributions from degree 37 to 360 
from the topographic/isostatic model (see Rapp and Pavlis, ibid. Section 3.4). Tn the new case we 
decided to reduce the maximum degree of the T2 model to 9 and then use the topographic/isostatic 
model from degree 10 to 360. This procedure makes our adjustment procedure somewhat more 
correct since we do neglect the correlation between the fill-in anomalies and the GEM-T2 
coefficients in our adjustment process. 

The second modification, then, took place by first generating a set of 30' anomalies, using 
program F431 using the GEM-T2 potential coefficients to degree 9 and the potential coefficients 
implied by the topographic/isostatic model from degree 10 to 360. These anomalies replaced the 
“fill-in” values after the China data had been merged. The code assigned to these fill-ins was 
3003. The standard deviations used were identical in the new case and the older case. This data 
was merged with the anomalies derived from satellite altimeter data using the same process 
described in Rapp and Pavlis (ibid, Section 3.5). The location of the 30' anomalies, with an 
indication of their origin is given in Figure 24. This figure can be compared to Figure 4 of Rapp 
and Pavlis where the only change will be seen in China. Statistics on the number of anomalies that 
are in the 30' data file (TS0040.DG30X30.MRGD.GEMT209.TI10360.CHINA) are given in 
Table 19. 
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Table 19. Statistical Information on the Various 30' Anomalies in the Merged Data File. 
Anomaly Units are mgal. 



30’ 

Terrestrial 

Altimeter 

Derived 

1° 

Split Up 

Fill-Ins 

Number of values 

45,166 

136,270 

30,187 

47,577 

Percentage by area 

17.51 

64.81 

7.95 

9.74 

Minimum value 

-196 

-284 

-175 

-243 

Maximum value 

391 

309 

214 

310 

Mean value 

3.8 

- 1.7 

6.2 

-0.3 

RMS value 

33.8 

25.0 

40.1 

27.3 

RMS standard dev. 

8.5 

3.5 

33.4 

36.0 


Using the GEM-T2 model, its error covariance matrix, and the revised 30' data file the 
combination solution was carried out using the same procedures as described in Rapp and Pavlis 
(ibid). The solution could be called a OSU89B type with a revised 30' anomaly data set that was 
consistent with the 1° normal equation, development described in Section 2. The 30' anomaly 
weighting procedure used in this analysis was identical to that used for OSU89B; that is the 
original standard deviation was multiplied by 2 and the resultant value restricted to fall in the range 
14-27 mgals. In Table 20 statistical information for this combination solution is given. 

Table 20. Statistical Information on the OSU89B Type Combination Solution with the OSU9 1 
30' Data File. 


Quantity 

Value 

Standard Deviation Range 

as 

Number of v^g > 7 mgal 
k 

RMS v, mgal 
Maximum Ivl 

14-27 mgal 

4.51 

14,492 

1.096 

3.21 

25 


The location of 14492 30' anomaly residuals exceeding 7 mgal in absolute value is shown 
in Figure 25. This number is substantially higher than the 4130 values shown in Fig. 5 of Rapp 
and Pavlis (ibid). The primary reason for this is the use, in our current solution of the fill-in 
anomalies from GEM-T2 to 9 while the OSU89B solution used GEM-T2 to degree 36. The 
patterns in Figure 25 are primarily those in regions where the fill in anomalies are used. 

After the least squares adjustment was completed a set of adjusted 30' mean anomalies was 
calculated. These anomalies were then converted to potential coefficients using the 
orthogonalization process described in Rapp and Pavlis (1990, Section 2.1). This led to a set of 
potential coefficients complete to degree 360. The coefficients from degree 51 to 360 will be used 
to augment the 2 to 50 coefficients found from the adjustment described in Section 4.0. However 
it still is of interest to compare the coefficients of the two solutions, up to degree 50, just to obtain 
a feeling for the differences to be seen from the two adjustment methods using two very different 
ways in which the altimeter data are treated. These comparisons are shown in Table 21 . From the 
table we see that the undulation difference of the two models is 48 cm. It will turn out that the 
accuracy of the OSU91A model is 25 cm to degree 50 (see Table 22) and the accuracy of the 
OSU89B type adjustment, to degree 50, is approximately 34 cm, so that the difference noted in 
Table 20 are consistent with the accuracy estimates of the two models. 
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Figure 25. Location of 14,492 30' x 30' cells where the anomaly residuals exceeded 7 mgal in absolute value. 





Table 21. Comparison of Two Potential Coefficient Models, To Degree 50, That Differ in 
Adjustment Procedure and Treatment of Altimeter Data 



At this point we have the set of potential coefficients to degree 50 from the adjustment 
described in Section 4.0 and a set of coefficients from degree 2 to 360 using the procedure 
described in the previous sections. In the next section we turn to the merger of the two models. 


63 













7.0 The Merger of the Two Potential Coefficients Sets Leading to the OSU91A 

Potential Coefficient Model 

We now simply take the potential coefficients to degree 50 from the combination solution 
of Section 4.0 and add the coefficients from degree 51 to 360 of the solution described in Section 
6. This merger leads to the OSU91A potential coefficient model. 

The standard deviations of the coefficients are also established in two steps. For the 
coefficients to degree 50 the standard deviations are based on the square root of the diagonal 
elements of the FYS 10.W96GW solution. The standard deviation of the coefficients from degree 
51 to 360 are exactly the same as the OSU89B coefficients because the anomaly accuracy estimates 
remain unchanged. The accuracy estimates for these higher degree coefficients are based on a 
propagated error and a sampling error as described in Section 2.3 of Rapp and Pavlis (ibid). Tabic 
22 shows the geoid accuracy, on a sphere of radius a, for the OSU89B solution and from the 
OSU91A solution. 

Table 22. Geoid Undulation Commission Error, by Spherical Harmonic Degree, for OSU89B and 
OSU91A. Units are cm. 


[ OSU89B 

OSU91A 

By Degree 

Cumulatively 

By Degree 

Cumulatively 

2 

0.2 

0.2 

0.2 

0.2 

6 

1.8 ! 

2.8 

1.3 

2.2 

10 

3.9 

7.6 

2.4 

5.0 

20 

5.4 

17.6 

3.6 

10.6 

30 

5.5 

24.8 

4.3 

16.8 

50 

4.6 

33.5 

3.0 

24.8 

75 

3.7 

39.3 

3.7 

32.3 

100 

3.2 

42.9 

3.2 

36.5 

180 

2.2 

48.7 

2.2 

43.2 

360 

1.3 

53.6 

1.3 

48.7 


Up to degree 50 the OSU91A has slightly smaller standard deviations than OSU89B. 
Beyond degree 50 the two solutions have identical accuracies at each degree. Overall, OSU9 1 A, to 
degree 360 has a commission error, only, of 48.7 cm vs 53.6 cm for OSU89B. 

The standard deviation of the anomaly degree variances, on a sphere of radius a, are plotted 
in Figure 26. Also shown on the plot are the anomaly degree variances (c n ) and the c n values 
implied by the Kaula rule for the decay of the fully normalized potential coefficients (Rapp and 
Pavlis (ibid, eq. 74). The anomaly power spectrum of OSU89B and OSU91A are very similar as 
can be seen by comparing Figure 26 with Figure 12 of Rapp and Pavlis (ibid). From Figure 26 we 
see that the signal to noise ratio becomes 1 near degree 260. Rapp and Pavlis (1990) have argued 
that coefficients above the degree at which the ratio becomes one should be retained in the solution. 
This also speaks to setting the error degree variances above degree 260 to be equal to the signal 
degree variances since the latter are smaller. In this case the commission error, now on a sphete 
whose mean radius is 6371 km, is 51.9 cm, which, combined with the omission error (based on 
the Tscheming/Rapp model with the Jekeli parameters) of 24.1 cm, yields a total point geoid 
undulation error of 57.2 cm. Recall this is a global average and can be poorer or better depending 
on the data availability in the area. One also should note in Figure 26 the break in the errors at 
degree 50. This break is caused by the inconsistency in the error estimates from the two 
combination models. We considered several techniques for artificially avoiding the break and 
finally decided to leave the errors as they were from the solutions. 

In the next section we evaluate the OSU91A model, and several preliminary versions, by 
comparisons with a variety of data types. 
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Figure 26. Anomaly signal and error degree variances impliecTby OSU91A and Kaula's rule. Values refer to the sphere of 
radius 6378137 m. 






8.0 Model Comparisons and Evaluations 

In this section we will evaluate some of the preliminary models, as well as the final 
geopotential model, OSU91A, developed for this report. In addition, for some evaluations, we 
will incorporate other recently published or developed geopotential models. 

8.1 Orbit Accuracy Assessment 

In this section we discuss orbit fits or orbit residuals when a geopotential model is used in 
an orbit determination process. In principle the smaller the residual fit the better the model. 
However the situation is not as simple as this because a potential coefficient model is only one part 
of a package used in the orbit estimation process. For example, one needs a set of station 
coordinates that should be consistently estimated with the geopotential model- In the tests to be 
described here the station coordinates, tidal models, etc., are held fixed at previously determined 
values. In some cases drag, radiation, etc. parameters may be estimated. Rapp and Pavlis (1990, 
Table 10) report, orbit accuracy assessments for a variety of models being evaluated at that time. 

The first group of tests were carried out by Fell (1991, private communication) using 
Doppler data acquired from 45 tracking station for Geosat and NOVA 3 satellites. The results of 
the tests are given through "station navigation statistics which are the weighted RMS in the radial 
and tangential directions computed from adjusting the geodetic position of each station on a single 
pass basis to best fit the estimated ephemeris for the satellite which is held fixed in the adjustment" 
(Sloop, 1991, private communication). 

These statistics are given for several models in Table 23. 

Table 23. Navigation Statistics for Various Geopotential Models Based on Doppler Tracking of 
Geosat Units are meters. 



RMS 1 

Residuals 

Model 

Radial 

Tangential 

GEM-T1 

2.2 

2.3 

GEM-T2 

2.3 

2.2 

OSU89A 

2.4 

2.2 

OSU89B 

2.4 

2.2 

S10.W24GW 

2.4 

2.2 

S10.W96GW(91A) 

2.4 

2.2 


This table does not indicate much sensitivity to the geopotential models tested. In Table 24 
statistics for the NOVA 3 satellite are given. 

Table 24. Navigation Statistics for Various Geopotential Models Based on Doppler Tracking of 
NOVA 3. Units are meters. 



RMS Residuals 

Model 

Radial 

Tangential 

GEM-Tl 

1.8 

2.5 

GEM-T2 

1.4 

1.6 

OSU89A 

1.3 

1.6 

OSU89B 

1.3 

1.6 

S10.W24GW 

1.3 

1.3 

S10.W96GW(91A) 

1.3 

1.4 1 
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Again the results, except for GEM-T1, are all about the same indicating that the preference between 
the models can not be determined from the table. On the other hand the combination solution that 
are represented here do not show a deterioration in fit over a satellite alone model. This, in itself is 
encouraging. 

Additional orbit fits were carried out at the Goddard Space Flight Center by Klosko f 1 99 1 , 
private communication) with the 1/24 and 1/96 weighting procedures. Results using "frozen’’ data 
sets are given in Table 25. 

Table 25. RMS of Fit for Selected Geopotential Models On Selected Satellite Using Laser 
Tracking and Frozen Data Sets. Units are cm. 



Model I 

Satellite 

GEM-T2 

W24GW 

W96GW 

Lageos 

5.2 

5.4 


BEC 

28.5 

26.5 


Geos- 1 

20.9 

22.1 

21.7 

Geos-2 

50.9 

48.3 

46.9 

Geos- 3 

17.9 

19.6 

19.2 

Ajisai 

10.9 

11.7 

11.6 

Starlette 

14.7 

17.0 

17.0 


These tables seem to give a consistent picture that indicates the W96GW (OSU91 A) model 
gives slightly better fits than the W24GW solution. The results are mixed when compared against 
GEM-T2. In some cases the T2 models gives better results while in others the OSU models are 
slightly better. 


Putney et al. (1991) have reported data fits in orbit tests with independent data with several 
geopotential models. The TEG-2 model was developed at the University of Texas at Austin 
(Tapley, private communication, 1991). Results for laser observation on three satellites are given 
in Table 26. 

Table 26. RMS of Fit for Selected Geopotential Models On Three Laser Tracked Satellites (Putney 
et al., 1991). Units are cm. 


Model 

IHSI&239H 


Starlette 

GEM-T2 

3.05 


13.28 

GEM-T3 

2.91 

9.99 

11.77 

TEG-2 

2.93 

9.77 

12.04 

OSU91A 

2.97 


13.17 


These results show that OSU91A gives fits similar to the other models tested but it is always 
slightly poorer (2% to 12%) than the GEM-T3 model. 

Putney et al. (ibid) also described orbit fit tests with SPOT-2 Doris data. Selected results 
from Putney et al. are shown in Table 27. 
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Table 27. RMS Fits To Spot-2 Doris Data. Units are mm/sec. 



Arc 1 

Model 

503 

702 

GEM-T2 

6.1 

6.6 

GEM-T3 

3.5 

3.4 

! TEG-2 

4.0 

4.2 

OSU91A 

4.7 

3.7 


This table shows that the OSU91A has made a significant improvement over the GEM-T2 
model for both arcs. The fit for arc 503, is significantly better and marginally better for arc 702 
with the GEM-T3 model as opposed to OSU91A. The results indicate that the OSU combination 
model has yielded better results than the initial model but other models are better in these 
comparisons. 

From the tests reported in this section we conclude that OSU91A performs well in the orbit 
fit tests but it is not the best model for the satellite orbit fits. Based on the results of Putney et al. 
(ibid) GEM-T3 gives better fits by 2 to 25%. Whether this is important or not must be judged in 
the context of the other tests to be reported in the next sections. 

8.2 Model Undulation Comparisons at Doppler Positioned Stations 

The principle in this comparison is to calculate the geoid undulation from the geopotential 
model and compare the value to that implied by an ellipsoid height from Doppler positioning, and 
an orthometric height of the station. Before the comparison can be done the Doppler heights must 
be put into a geocentric reference system that is properly scaled. The procedures followed for this 
report are identical to those described by Rapp and Pavlis (1990, Section 4.2). The results are 
given in Table 28 where the mean difference (Doppler minus model), the standard deviation of the 
difference, and the number of accepted stations are given. Two models, GEM-T3 and TEG-2, 
were augmented in these comparisons, with the coefficients from degree 51 to 360 from the 
OSU91A model. 

Table 28. Comparison of Geoid Undulations from Doppler Positioning with Geoid Undulations 
from Potential Coefficient Models 



Mean 

Standard 

Number of 

Model 

Difference, m 

Deviation, m 

Stations 

OSU89B 

.14 

±1.60 

1782 

OSU91A 

.15 

1.58 

1802 

GEM-T3 

.12 

1.57 

1788 

TEG-2 

.10 

1.77 

1767 


These results indicate a slight improvement of OSU91A over OSU89B. The T3 comparisons 
shows 14 less stations (than 91A) used in the comparisons while the TEG-2 model has a 
significantly larger standard deviation. 

One can also test to see if the undulation discrepancies in these comparisons have a 
correlation with elevation. This is done by making a linear fit to the residuals and examining the 
height dependent term. Table 29 shows results of these comparisons including the RMS residual 
(after adjustment). 
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Table 29. Bias and Slope Fits to Undulation Residuals as a Function of Elevation for the Global 
Station Set 


Model 

Bias, m 

Slope, m/km 

RMS 

Residual, m 

OSU89B 

.059 ± .104 

-.008 ± M6 

.226 

OSU91A 

.034 ± .099 

.037 ± .082 

.216 

GEM-T3 

.006 ± .105 

.038 ± .087 

.230 

TEG-2 

-.131 ± .140 

.159 ± .116 

.305 


This table shows that the 91 A gives smaller residuals than the 89B model. The slope is negligible 
for 89B, 91 A, and T3 and at the edge of being significant for TEG-2. The RMS residual is 41% 
higher than that found for 91 A. 

We conclude by noting the 91 A model and the augmented T3 model appear to have 
comparable accuracy although 91A would have the edge because more stations are accepted for 
use, and the RMS residual, from Table 29, is smaller. 

Recently Shibuya (1991, private communication) made available the ellipsoidal height of a 
station (the S point) at Breid Bay, East Antarctica whose position was determined in the WGS84 
terrestrial reference frame using the precise DMA ephemeris. Also available (Shibuya, Fukuda, 
Michida, 1991) was a mean sea level height for the station. Taking into account the significant 
(-2.0 m) sea surface topography in this area, and referencing the undulation to the equatorial radius 
of 6378136.3 m, the geoid undulation was estimated as 22.8 m. Using OSU91A the geoid 
undulation was 22.2 m, an excellent agreement. A nearby Geosat track, again after correction for 
sea surface topography, implied a geoid undulation of 21.6 m. These two values are in excellent 
agreement. 

8.3 Comparison with Undulations and Undulation Differences at 
GPS/Orthometric Height Stations 

In the past several years GPS measurements have been made to accurately position 
stations, primarily in a relative sense. The positions are generally determined by fixing one station 
through coordinates defined in a terrestrial reference frame. This frame might be in a geodetic 
datum (e.g. NAD83) or in a space related frame defined by satellite laser ranging or VLBI 
measurements. Rapp and Pavlis (1990, Section 4) describe four GPS/leveling traverses that will 
also be analyzed here. In addition Rapp and Kadir (1988) describe an area network of GPS 
stations in the State of Tennessee that will also be examined here. The Tennessee network 
consisted of 49 stations distributed uniformly across the State of Tennessee. The positions were 
not placed in a geocentric system so that systematic undulation differences are to be expected 
between the geoid undulations implied by the potential coefficient model and the GPS ellipsoidal 
height/leveling orthometric height. A complete description of other traverses may be found in 
Rapp and Pavlis (1990, Section 4.3). 

Table 30 shows the mean differences (model minus GPS/leveling) between the two 
undulation estimates. 
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Table 30. Mean Undulation Difference at GPS Stations. Units are cm. 




' ”Tj 

fraverse or Networks 



Europe 

Canada 

Australia 

Scandinavia 

Tennessee 

Num of Stations 

60 

63 

38 

46 

49 

Model 

OSU89B 

26 

-38 

-3 

29 

149 

OSU91A 

40 

4 

-49 

33 

157 

GEM-T3 

35 

-16 

-18 

28 

183 

TEG-2 

78 

43 

61 

64 

262 


This table shows variations between the models and the traverse/network being tested. The largest 
mean difference occurs for Tennessee which has not been placed in a geocentric system. 

Table 31 shows the standard deviation (mean removed) of the undulation difference for 
the same traverse/network shown in the previous table. 

Table 31. Standard Deviation of the Undulation Difference at GPS Stations in Five Areas. Units 
are cm. 



1 

fraverse or Networks 

Model 

Europe 

Canada 

Australia 

Scandinavia 

Tennessee 

OSU89B 1 

32 

39 

34 

32 

30 

OSU91A 

33 

36 

35 

32 

21 

GEM-T3 

42 

36 

35 

48 

23 

TEG-2 

70 

57 

68 

75 

73 


Table 31 shows that 91A is slightly better than 89B although a significant improvement takes place 
for the Tennessee area. GEM-T3 is comparable to the other OSU solutions in Canada, Australia, 
and Tennessee and poorer for the European and Scandinavia traverse. The TEG-2 model performs 
poorer than the other models in these tests. 

The next comparison is a relative comparison where undulation differences are compared. 
Such a comparison removes or reduces long wavelength error in the data. In these computations 
the differences are computed between adjacent stations in the traverse. In the Tennessee networks 
comparisons are made provided the station spacing is less than or equal to 55 km. This distance 
limits the comparisons to adjacent stations in the network. The results are given in terms of RMS 
differences and in ppm (parts per million) by dividing the RMS difference by the length of the line. 
The values are given in Table 32. 

Table 32. Relative Geoid Undulation Comparisons for GPS/leveling and Model Results 



Traverse or Networks 

Model 

Europe 

Canada 

Australia 

Scandinavia 

Tennessee 


RMS 

ppm 

RMS 

ppm 

RMS 

ppm 

RMS 

ppm 

RMS 

ppm 

OSU89B 

25 

3.5 

9 

6.6 

22 

5.1 

24 

4.3 

28 

4.9 

OSU91A 

23 

3.6 

9 

6.6 

22 

5.3 

25 

4.4 

26 

3.9 

GEM-T3 

24 

3.7 

9 

6.6 

22 

5.3 

26 

4.6 

26 

3.9 

TEG-2 

30 

5.2 

10 

6.8 

26 

5.6 

33 

5.8 

57 

6.0 


The results in Table 32 indicate the first three models listed are of comparable accuracy while the 
TEG-2 model performs somewhat poorer. 
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In judging these differences it is of interest to compare the results with what is expected 
from the errors given for the coefficients of the OSU91A model. To do this we write the 
commission error for a undulation difference computed from potential coefficients to degree N and 
perfect gravity data in a spherical cap of radius to Yo as follows (Christodoulidis, 1976, eq. (135)): 

C?Q = 4X Q£(Yo)EnS n+2 [l - Pn(cOS\|/pQ)] 

2fn=2 (8.1) 

where: R is a mean Earth radius: 

y is an average value of gravity 

fo, is the Molodensky truncation coefficient at degree n; 

£n is the error anomaly degree variances; 
s is a factor to refer results to the mean sphere. 

The omission or truncation error would be written as: 

oo 

= Q2(Yo)CnS n+2 [ 1 - P^cosypq)] 

2 r N+1 (8.2) 

where c n are the anomaly degree variances implied by a model. 

If only potential coefficient information is to be used, Yo = 0 and 

Qn<¥o = 0) = S A T (83) 


Equations (8.1) and (8.2) now become: 


hr e 2 s n+ 2 [l . P n (coSYPQ)] 

T n= 2 (n -l) 2 (8.4) 


t pq = 2 T-X 7 bz n+2 t 1 • p n(cosYPQ)] 

T n +1 (n - lr (8.5) 

Letting R = 6371000 m, y = GM/R 2 = 9.820 ms* 2 , and using the Jekeli parameters in the 
Tscherning/Rapp degree variance model equation (8.4) and (8.5) have been evaluated. The £n 
values have been based on the actual error degree variances of OSU91 A to degree 260, and the 
actual anomaly degree variances of 91 A from degree 261 to 360. The individual results and total 
results are given in Table 33 for the OSU91 model. In addition the total error has been expressed 
in parts per million of the distance. 
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Table 33. Commission and Omission Error in the Computation of Geoid Undulation Differences 
from the OSU91A Potential Coefficient Model 


Linear 

Separation 

Angular 

Separation 

HEBHH 

Truncation 

Error 

Total 

Error 

H 


(TOO 

0 cm 

0 cm 

0 cm 

0 cm 


.09 

9.2 

10.3 

13.8 

13.8 


.18 

18.1 

17.9 

25.5 

12.7 


.27 

26.6 

22.7 

35.0 

11.7 

40 

.36 

34.4 

25.3 

42.7 

10.7 

50 

.45 

41.3 

26.2 

48.9 

9.8 

70 

.63 

52.4 

24.8 

57.9 

8.3 

90 

.81 

59.6 

22.6 

63.7 

7.1 

100 

.90 

62.0 

22.1 

65.8 

6.6 

200 

1.80 

71.8 

23.1 

75.4 

3.8 

400 

3.60 

77.1 

23.6 

80.6 

2.0 

600 

5.40 

77.7 

23.6 

81.2 

1.4 

800 

7.19 

77.7 

23.5 

81.2 

1.0 

1000 

8.99 

77.8 

23.5 

81.2 

0.8 

1600 

14.39 

77.4 

23.5 

80.9 

0.5 

2000 

17.99 

77.4 

23.5 

80.9 

0.4 

10000 

89.93 

77.4 

23.5 

80.9 

.08 


The results from Table 33 can be compared to results from actual comparisons as represented in 
Table 32. Note that the errors given in Table 32 are errors in all data used in the comparison: 
geopotential model; GPS; leveling. Even then the actual results are much better than is expected 
from the projected OSU91 model errors shown in Table 33 noting the average line length in Table 
32 is 30 km. 

8.4 Comparison with Undulation and Undulation Differences at Two VLBI Sites. 

Forsberg and Madsen (1990) have pointed out that precise ellipsoidal and orthometric 
heights exist at the VLBI sites which have been tied into the geocentric terrestrial reference system 
FSC-2. Knowing the ellipsoidal height in the geocentric system and the orthometric height a 
"standard" geoid undulation can be defined. The coordinates of the two stations are as follows: 

Tromso 69°40’N, 18°56'E 

Onsala 57°24’N, H°56'E 

The standard undulation, undulation from a variety of models, and undulation differences relative 
to the standard are given in Table 34 for Tromso and Table 35 for Onsala. 

Table 34. Undulation and Undulation Discrepancies at the Tromso VLBI Station 


Solution 


Difference w.r.t.VLBI 

VLBI 

31.11 m 

— 

NKG89* 

31.45 

34 cm 

OSU89B 

31.67 

56 

OSU91A 

31.22 

11 

GEM-T3 

31.48 

37 

TEG-2 

31.36 

25 


* detailed geoid solution; see Forsberg and Madsen. 
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Table 35. Undulation and Undulation Discrepancies at the Onsala VLBI Site 


Solution 

Undulation 

Difference w.r.t.VLBI 

VLBI 

36.54 m 

— 

NKG89* 

36.28 

-26 cm 

OSU89B 

36.65 

11 

OSU91A 

36.71 

17 

GEM-T3 

36.54 

0 

TEG-2 

36.32 

-22 


These two tables show different levels of agreement, all of which are fairly good. It is interesting 
to note that the OUS91A and TEG-2 model both show better agreement with the "standard" than 
the detailed geoid, NKG89. 

A final comparison can be made by calculating the undulation difference, Tromso - Onsala, 
from VLBI and from the models. These results are given in Table 36. 

Table 36. Geoid Undulation Difference Comparisons Between the Tromso and Onsala VLBI Sites 



Value 

Difference w.r.t.VLBI 

VLBI 

-543 cm 

— 

NKG89 

-483 

60 cm 

OSU89B 

-498 

45 

OSU91A 

-549 

6 

GEM-T3 

-506 

37 

TEG-2 

-496 

47 


The results from this table show a great variety of differences, with OSU91A clearly giving the 
best agreement with the standard. 

8.5 Comparison with Geoid Undulations Derived from Geosat Altimeter Data 

In this section we discuss the geoid undulations derived from Geosat altimeter data with the 
geoid undulations derived from the various geopotential models. Such comparisons were 
described in Rapp and Pavlis (1990, Section 4.4). The basic principle in finding the geoid 
undulation from satellite altimeter data is to assume that sea surface topography (Q is known, and 
to remove it from a sea surface height that has been determined from an accurate satellite 
ephemeris. Specifically we have: 

N = SSH - C (8-6) 

For the calculations to be described here the sea surface heights have been calculated from the 
orbits developed as part of the adjustment process described in Section 2 and Section 4. The sea 
surface topography was the (10, 10) model, based on the 1/96 altimeter weighting, as described in 
Section 5. In the comparisons the Geosat data consisted of 768286 points in the following time 
period: 870220 (11 hr 8 min), 870307 (23 hr 55 min). This 17 day time frame corresponds, 
approximately, to the Exact Repeat Mission 7. This time period has been used because it 
corresponds to the summer time in the Antarctic region so that data, in normally ice covered areas, 
is available. 

In carrying out these comparisons it is important to recognize the role of the permanent 
Earth tides (Rapp, 1989, Rapp et al., 1991). In our case, the Geosat surface is in a "mean" 
system. The geoid undulations (N n ) that are computed from the potential coefficient models are 
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placed in the non-tidal system which is consistent with the system used for the development of 
GEM-T2. We then define the undulation discrepancy as follows: 

d = N n - 0.257 (| sin 2 <l>-l) - (SSH - C) (8J) 

where the second term on the right hand side of (8.7) follows from eq. (18) of Rapp (1989) with 
k 2 = 0.3. Table 37 gives statistical information on "d" for several geopotential models. This table 
is analogous to Table 8 of Rapp and Pavlis (1990). For consistency reasons the sea surface 
topography was kept at the OSU91 degree 10 model. However one test was carried out with a 
degree 10 (TEG-2) model supplied by Shum (1991, private communication). The geopotential 
model is defined by the model identified augmented by the OSU91A model from degree 5 1 to 360 
except in the case of OSU89B which is itself complete to 360. 

Table 37. Results from the Comparison of Model and Geosat (ERM 7) Implied Geoid 
Undulations. Units are m. 


Differences 

OSU89B 

OSU91A 

GEM-T3 j 

TEG-2 

TEG-2* 

Minimum 

^34 

-6.92 

-6.52 

-8.67 

-7.71 

Maximum 

5.53 

4.54 

4.49 

10.03 

10.86 

Mean 

0.005 

0.008 

0.019 

0.033 

0.057 

Std. Dev. 

0.529 

0.342 

0.485 

0.739 

0.785 

No > 1.5 m 

14291 

4505 

16204 

43181 

48298 

No > 2.0 m 

4008 

1456 

5176 

22280 

25929 

No > 4.0 m 

82 

39 

128 

1774 

3196 


* with TEG-2 SST model. 

From this table we see that the 89B model gives better results (53 cm vs 61 cm) than presented in 
Table 8 of Rapp and Pavlis (ibid). This is presumably due to the improved Geosat orbits described 
in this paper. The OSU91A model shows a significant improvement oyer 89B (± 34 cm vs 53 cm) 
with a considerable reduction (14291 to 4505) in the number of discrepancies > 1.5 m. The 
augmented GEM-T3 model performs better than the TEG-2 model in these tests but still does not 
give as good a fit as the OSU91A model. The last two columns of Table 37 differ in the sense that 
the last column uses the TEG-2 SST model while the second from the last uses the OSU91 SST 
model. The comparisons are better when the OSU91 model is used. This may be due to the error 
in the TEG-2 SST model in geographic areas for which data was not included in the solution. In 
assessing these comparisons one must recognize that the orbits being used are those consistent 
with the OSU91A gravity field up to degree 50. Since the orbit is dependent on the gravity field, it 
would be appropriate to carry out these comparisons with Geosat orbits that are consistent with the 
GEM-T3 or TEG-2 model. No such tests were carried out because such orbits were not available. 

It is appropriate to point out that the standard deviation of 34 cm associated with the 
OSU91 model is smaller than would be expected from the discussion near Table 22 where it was 
argued that the total undulation error, at a mean radius of 6371 km, is 57 cm, considerably more 
than the 34 cm found from the Geosat comparisons. However the 57 cm is a global estimate and 
the results for the ocean area would be slightly smaller but not enough to account for the apparent 
pessimistic accuracy estimate obtained through formal error propagation. 

We next examine the location of the residuals that exceeded 1.5 m for the 91 A (Figure 27), 
GEM-T3 (Figure 28), and TEG-2 (Figure 29) models. A similar figure for the OUS89B model is 
Figure 8 in Rapp and Pavlis (ibid). 
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Figure 27, for the OSU91A model shows that the discrepancies occur primarily in areas of 
high frequency where contributions from the harmonics above degree 360 may be significant. 
Figure 28, for the augmented GEM-T3 model shows the high frequency discrepancies noted for 
the OSU91A model, but also discrepancies in the Mediterranean and Caspian Seas as well as large 
patches of discrepancies below -60°. Figure 29 for the augmented TEG-2 model (with the OSU91 
SST model) shows the high frequency effect as well as substantial discrepancies in the 
Mediterranean and Caspian Seas and large patches of discrepancies below -60° latitude. In addition 
track related discrepancies are obvious north of Australia and west of South America. 

The large patterns of residuals for the GEM-T3 and the TEG-2 model is apparently 
traceable to editing criteria that deleted data from these areas in the analysis leading to GEM-T3 and 
TEG-2. This criteria specifically deleted data below -60° latitude, in shallow water areas, and in 
areas where the tide model was not defined on the GDR. These editing criteria would seem to be 
too 
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Figure 27. Location of 4505 Points on Geosat ERM7 Where the Corrected Geoid Model (OSU91 A) Minus the Corrected Sea 
Surface Height Is > 1.5 m. 
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Figure 29. Location of 43181 Points on Geosat ERM7 Where the Corrected and Augmented Geoid Model (TEG-2) Minus 
the Corrected Sea Surface Height Is > 1.5 m. 




restrictive and have led to a gravity model that does not fit the data in areas where the data was not 
used in the solution. No surprise. 

Another set of comparisons were made when data below -60° latitude, and in the 
Mediterranean area (30° < <]) < 46°, 0° < X, < 38°) were excluded from the comparisons. These 
results are shown in Table 38 where 668534 points were compared. 

Table 38. Altimeter/Model Geoid Undulation Comparisons Deleting Data Below -60° and in the 
Mediterranean Area. Units are m. 



Model | 


OSU91A 

GEM-T3 

TEG-2 

TEG-2* 

Std. Dev. 
No £: 1.5 m 

0.34 

3656 

0.40 

6054 

0.53 

15672 

0.53 

14589 


* with TEG-2 SST model. 

From this table we see no change in the standard deviation for the OSU91A model as compared to 
the complete case given in Table 37. The results for GEM-T3 are improved significantly although 
the TEG-2 model shows residuals 2.5 times greater in number than GEM-T3. Additional tests are 
described by Rapp, Wang, and Pavlis (1991a) in which sea surface heights along Geosat tracks on 
the Mediterranean Sea were compared with the corresponding undulation from selected models. 
With the TEG-2 model the discrepancies reached 6.2 m while other models showed better 
agreement. Statistics on the agreement along two Geosat arcs are given on Table 39. 

Table 39. Standard Deviation of Geoid Undulation Differences for Two Geosat Tracks (ERM 7) 
in the Mediterranean Sea 


Model 

Track 90874 

Track 90888 

PGS4233 

45 cm 

242 cm 

TEG-2 

139 

625 

OSU91A 

20 

54 


Track 90874 starts south of Italy and goes to the African coast. Track 90888 crosses the Agean 
Sea and the Mediterranean. It is in an area of extreme variations of the geoid accounting for the 
poorer fit for the track as compared with track 90874. 

8.6 Conclusion 

This lengthy section has been developed to report comparisons of a number of geopotential 
models with selected data for evaluation purposes. The evaluations range from orbit fit tests to the 
comparison of undulation differences from VLBI/leveling data. The orbit fits with OSU91A 
showed reasonable results but for some tests (e.g. with DORIS data) not quite as good as GEM- 
T3. Many other orbit tests could be run and the results would depend on how the satellite data was 
weighted in the solution. The tests with the Doppler undulation data showed the 91 A and T3 
solutions are comparable and TEG-2 is poorer. The undulation comparisons at GPS/leveling 
stations showed 91A gives somewhat better results than GEM-T3 while TEG-2 is significantly 
poorer than 91A or T3. Undulation and undulation difference comparisons at two VLBI sites in 
Europe indicate good agreement with 91 A. 

One of the most revealing tests was with the altimeter implied geoid undulations. Using a 
complete Geosat ERM statistics on the undulation differences for several models were computed. 
The OSU91A model gave the best agreement which is to be expected to some extent since there 
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was Geosat data used in the solution to degree 50. Poorer comparisons were obtained with GEM- 
73 and TEG-2 although GEM-T3 was significantly better than the TEG-2 model. Examination of 
the location of larger residuals indicated areas for which data was edited out in the T3 and TEG-2 
data analysis. These results indicate careful consideration must be given to data editing since 
models will not fit well in areas where no data has been included in a solution. 
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9.0 Undulation Differences Between Models 

Rapp, Wang, and Pavlis (1991b, 1991c) have compared several geopotential models 
primarily in terms of global undulation differences. Selected results from these papers arc 
described in this section. 

We first compute the geoid undulation on a 1° x 1° grid from the potential coefficient model 
to degree 50. Table 40 gives the standard deviation of the global differences, Table 41 gives the 
differences for land (positive elevation) and Table 42 gives the differences for ocean (negative 
elevations) areas. 

Table 40. Standard Deviation of Geoid Undulation Differences on a Global Scale. Units are cm. 


Model 

91 A 

T3 


OSU91A 

GEM-T3 

47 

47 

106 

83 


Table 41. Standard Deviation of Geoid Undulation Differences in Land Areas. Units are cm. 


Model 

91A 

T3 

TEG-2 

OSU91A 

— 

73 

176 

GEM-T3 

73 

— 

131 


Table 42. Standard Deviation of Geoid Undulation Differences in Ocean Areas. Units are cm. 


Model 

91 A 

T3 

TEG-2 


— 

24 

65 


24 

— 

47 


The best agreement between the models takes place in the ocean areas. The ±73 cm difference 
(SD) between 91A and GEM-T3 in land areas can be attributed to different weighting procedures. 
The differences between TEG-2 and the OSU91A and GEM-T3 are due to different surface normal 
equations used, TEG-2 used a set of normal equations based on 1986 surface data while OSU91A 
and GEM-T3 used a set of normal equations based on 1990 data. Table 43 gives the magnitude of 
the largest discrepancies: 

Table 43. Maximum Absolute Geoid Undulation Difference Between Selected Geopotential 
Models to Degree 50. Units are cm. 


Model 

91A 

T3 

TEG-2 

OSU91A 

GEM-T3 

446 

446 

910 

623 


The 9.1 m difference between 91A and TEG-2 occurs in western South America near latitude -19°. 
An additional large difference of 9.0 m occurs in the eastern Mediterranean Sea near longitude 30°. 
The largest difference between the 91 A and GEM-T3 model occurs in the Himalaya Mountains. 
Figure 30 shows the undulation differences GEM-T3 minus OSU91A. The good agreement (± 24 
cm) in the ocean areas is clearly contrasted with the poorer agreement (±73 cm) in the land areas. 
Also clear from this figure is the poorer agreement below -69° and in the Mediterranean Sea. 

In Figure 31 we show the geoid undulation differences by degree for GEM-T3 and TEG-2 
with respect to the OSU91A model. TEG-2 shows significant differences after degree 10. At 
degree 10 there is a 5.8 cm difference between OSU91A and TEG-2 while the difference is 3.8 cm 
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LONGITUDE 



Figure 30. Geoid Undulation Differences, GEM-T3 Minus OSU91A, Band on a 3° x 3° Grid with Coefficients to Degree 50. 
Contour Interval is 40 cm. 








with GEM-T3. Table 44 gives the undulation differences by degree. Figure 32 shows the 
cumulative undulation difference up to the specified degree. 

Table 44. Geoid Undulation Differences By Degree With Respect to OSU91A. Units tire cm. 


Degree 

GEM-T3 

TEG-2 

2 

0.7 

2.4 

5 

2.9 

2.9 

10 

3.8 

5.8 

20 

6.7 

15.9 

30 

8.7 

17.9 

40 

8.5 

19.4 

50 

6.4 

18.4 

Cumulative 

47 

106 


The better agreement with GEM-T3 may be attributed, in part, to the OSU91A model development 
starting from the GEM-T2 model. 

This section has examined the geoid undulation differences between three geopotential 
models: OSU91A, GEM-T3, and TEG-2. All comparisons have been made to degree 50. We 
find that OSU91A agrees better with the GEM-T3 model than with TEG-2. The models agree best 
in the ocean areas with significant (9 m) differences in some areas. 

Similar comparisons could be made for gravity anomalies and for comparisons with the 
OSU89B potential coefficient model which is complete to degree 360 as is the OSU91 A model. 
Table 45 gives such comparisons. 

Table 45. Cumulative Undulation and Anomaly Differences Between OSU91A and OSU89B. 

Up to Degree 50 and 360. 



Maximum Degree 


50 

360 

Geoid Undulation 
Gravity Anomaly 

57 cm 
2.2 mgal 

59 cm 
3.5 mgal 


The small change in going from degree 50 to 360, in the geoid undulation, reflects the small 
change in gravity material used in the 91A solution from the 89B model. 
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10. Conclusions 

This report describes the computation of a geopotential model to degree 360, a sea surface 
topography model to degree 10/15, and adjusted Geosat orbits for the first year of the exact repeat 
mission(ERM). This study started from the GEM-T2 potential coefficient model and it's error 
covariance matrix and Geosat orbits (for 22 ERMs) computed by Haines et al. (1990) using the 
GEM-T2 model. 

The first step in the study followed the general procedures described by Denker and 
Rapp(1990). This procedure used a radial orbit error theory originally developed by 
Engelis(1987a). The Geosat data was processed to find corrections to the a priori geopotential 
model, corrections to a radial orbit error model for 76 Geosat arcs, and coefficients of a harmonic 
representation of the sea surface topography. The processing of Denker and Rapp(ibid.) was 
extended by the addition of surface gravity normal equations from a recent 1° x 1° gravity anomaly 
data set. Using this data strengthened the separation between the potential coefficients and the sea 
surface topography harmonic coefficients. The r esu lts of this analysis led to sea surface estimates 
that were more realistic than found in the Denker/Rapp study where the sea surface slope across the 
Pacific Ocean was unrealistic when judged against oceanographic data. The results presented in 
this report appear realistic and attest to the improvement made in the Geosat orbits based on the 
GEM-T2 model and consistent tracking station coordin ate s. The RMS residual and crossover 
discrepancies were about the same as before (20 cm and 18 cm respectively). In this analysis we 
included altimeter data in areas not used previously (e.g the Mediterranean Sea and areas of high 
frequency gravity signal). This led to improved gravity fields in these areas as well as a sea surface 
topography estimate in the Mediterranean area that seems realistic. Our final solution in this first 
stage of processing was selected after testing several different weighting options and testing 
(through the help of the Goddard Space Flight Center) the different solutions in orbit fit tests. The 
primary weighting considerations related to the altimeter data weights trying to take into account the 
duplicative gravity signature sensed by the 22 repeat orbits analyzed. The weights used for the 
surface gravity data were not studied in detail but were selected considering the original standard 
deviations of the original data and calibration factors computed in several test solutions. Error 
covariance matrices were computed for the potential coefficients as well as the sea surface 
harmonic coefficients. The potential coefficient errors were used to compute geoid undulation 
accuracies at specific degree as well as cumulatively. For example up to degree 10 the cumulative 
geoid undulation error was 5.0 cm which is 25% smaller than reported by Rapp and Pavlis lor the 
OSU89B potential coefficient model. 

Time variations of the sea surface topography were also studied through spherical harmonic 
expansions carried out at one month intervals. These variations were used to compute sea level 
changes on an ocean wide basis. In addition, specific areas in the Pacific Ocean were studied to see 
if the sea level changes reported by Cheney and Miller(1990) could be detected. Although the 
agreement was not perfect it was clear that the techniques developed for this study gave reasonable 
agreement with the Cheney and Miller results. Based on the monthly sea surface topography 
solutions, the spectral analysis of the time variations of the sea surface height were also studied. 
The analysis of sea surface topography in this report should be regarded as a first step with 
additional study needed. 

The second stage of geopotential analysis took place by carrying out a combination of the 
GEM-T2 coefficients with 30' gravity data derived from surface gravity data and anomalies 
obtained from altimeter data, primarily, GEOS-3 and SEASAT. This data was identical to that used 
by Rapp and Pavlis (ibid.). For areas lacking gravity information the gravity anomalies were 
derived from a potential coefficient model defined by the GEM-T2 model to degree 9 plus the 
coefficients from a topographic-isostatic model to degree 360. This was different from the 
approach used in the OSU89B model development where the GEM-T2 model was taken to degree 
36 plus the topographic information. The combination solution gave a set of coefficients complete 


86 


to degree 50 plus a set of adjusted 30' mean gravity anomalies. These anomalies were used to 
compute the potential coefficients to degree 360. The coefficients of this model to degree 50 were 
replaced by the coefficients to degree 50 found in the first stage of the analysis. By doing this we 
take advantage of the improved analysis of the Geosat data that was carried out starting from the 
GEM-T2 Geosat orbits. The merged set of potential coefficients were designated the OSU91A 
potential coefficient model. The standard deviations for the coefficients were computed from the 
error covariance matrix and the error propagation/sampling error considerations for the coefficients 
above degree 50. 

The 91 A potential coefficient model was evaluated by comparison with several external data 
types. Orbit fit tests were carried out at the Goddard Space Flight Center and the Naval Surface 
Weapons Center. These tests indicated that the 91 A model was comparable to other current 
potential coefficient models but other models could be better for some satellite applications. Tests 
were made comparing geoid undulations derived from space and levelling data and the 
corresponding quantities derived form the geopotential models. These tests indicated that the 9 1 A 
model was about comparable with the 89B model although in some cases (e.g. in the State of 
Tennessee and at two VLBI sites) better. The undulations computed from the 91 A model were also 
compared to the undulations implied by a complete ERM of Geosat data computed from the new 
Ohio State orbits. The standard deviation of the difference between the two undulation estimates 
was 34cm a considerable improvement over the OSU89B 53 cm standard deviation. The 34 cm fit 
was considerably better than fits obtained with the GEM-T3 or TEG-2 geopotential coefficient 
models (augmented by the OSU91A coefficients above degree and order 50). The most serious 
problem encountered with these other models is the large number of residuals in areas where 
altimeter data was edited out on the solutions. The specific areas include the Mediterranean Sea and 
the region below -60 latitude. Such data was included in the OSU91A model development. 

The analysis in this report has shown how we can determine a high degree spherical 
harmonic model combining the best aspects of two different analysis techniques. The selection of 
weights for the data combination were made to give an overall fit to different data types. In tins 
sense it is a compromise solution. Tests however indicate the compromise has not significantly 
degraded the solution for any of the data type being tested. Additional study is need to learn about 
optimum weighting procedures in the sense that data fits tests are needed with a variety of data. 

There are a number of things that could be done to improve the solutions described in this 
report. For example, we have used the wet tropospheric correction that was given on the original 
Geosat GDRs. Literature has shown that these estimates may be in error by an amount that is 
significant with respect to the accuracy we are trying to achieve. Study is therefore warranted on 
the influence of the wet tropospheric correction errors on the parameters estimated in this study. 
We are also concerned about the appropriate representation of sea surface topography. As pointed 
out in Denker and Rapp(ibid.) the harmonic coefficient representation may not be the most 
appropriate when dealing with a non-global data set. Hwang(1991) has recently explored an 
alternative approach to this representation which may be much more satisfactory than the current 
modeling procedure. This modeling needs to be tested. 

The proper reduction of the surface gravity material to the geoid seems to be an important 
factor in the analysis where cm accuracy in geoid computations is being attempted. Our studies 
have used an assumption about gravity and elevation correlations to compute certain corrections 
based on elevation data. We need to investigate this assumption. In addition we need to improve 
the elevation data models so that more reliable correction terms can be computed. Specifically this 
calls for the improvement of the ETOP5U 5'x5' elevation model that is distributed by the National 
Geophysical Data Center. 

This report has described the error analysis that has led to the accuracy estimates for all the 
coefficients to degree 360. However this analysis from degree 51 to 360 is subject to a number of 
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assumptions that need to be evaluated. We are concerned that the accuracy estimate projected for 
the geoid undulation is significantly higher than that found by comparison with altimeter data in sea 
and undulation differences on land. The development of the higher degree coefficients throughout 
the orthogonality relationships also rests on assumptions related to the development of optimum 
desmoothing factors. The solutions being used in this report date back to studies made in the early 
1980's. Newer studies are needed in this area. The most recent attempt at doing this is described 
by Sacerdote and Sanso(1991). 

Although the analyses reported in this paper have shown improved ways to estimate the 
Earth's gravitational potential significant work is need to improve the modelling effort. We need 
additional surface gravity data for many regions now lacking data. We need to strive to get a global 
land coverage of 30' mean gravity anomalies. We need an improved analysis of satellite altimeter 
data to determine mean anomalies in the ocean areas. This new solution would take into account the 
improved Geosat orbits we now have as well as the new reference geopotential developed as part 
of this study. And we need the analysis in several areas as noted in the paragraphs above. 
Combining all this information together can help us reach an improved determination of the Barth's 
gravity field and the ocean surface through the altimeter data. 


88 


REFERENCES 


Brigham, E.O., The Fast Fourier Transform and Its Applications, Prentice Hall, Englewood 
Cliffs, N.J., 1988. 

Cheney, R. and L. Miller, Recovery of the Sea Level Signal in the Western Tropical Pacific from 
Geosat Altimetry, J. Geophys. Res. 95, C3, 2977-2984, 1990. 

Christodoulidis, D., On the Realization of a 10 cm Relative Oceanic Geoid, Report No. 247, Dept, 
of Geodetic Science and Surveying, The Ohio State University, Columbus, 1976. 

Cizek, V., Discrete Fourier Transforms and Their Applications, ADAM Hilger Ltd., Bristol and 
Boston, 1986. 

Denker, H. and R. Rapp, Geodetic and Oceanographic Results from the Analysis of 1 Year of 
Geosat Data, J. Geophys. Res., 95, C8, 13,151-13,168, 1990. 

Denker, H., Radial Orbit Error Reduction and Sea Surface Topography Determination Using One 
Year of Geosat Altimeter Data, Report No. 404, Dept, of Geodetic Science and 
Surveying, The Ohio State University, Columbus, 69A, 1990. 

Despotakis, V., The Development of the June 1986 1° x 1° and the August 1986 30' x 30' 
Terrestrial Mean Free-Air Anomaly Data Bases, Internal Report, Dept, of Geodetic 
Science and Surveying, The Ohio State University, 1986. 

Elliott, D.F., and K.R. Rao, Fast Transforms, Algorithms, Analyses, Applications, Academic 
Press, N.Y., London, Paris, 1982. 

Engelis, T., Radial Orbit Error Reduction and Sea Surface Topography Determination Using 
Satellite Altimetry, Report No. 377, Dept, of Geodetic Science and Surveying, The 
Ohio State University, Columbus, 1987a. 

Engelis, T., Spherical Harmonic Expansion of the Levitus Sea Surface Topography, Report No. 

385, Dept, of Geodetic Science and Surveying, The Ohio State University, Columbus, 
1987b. 

Forsberg, R. and F. Madsen, High-Precision Geoid Heights for GPS Levelling, presentation at 
GPS-90 Symposium, Ottawa, September 1990. 

Haines, B.J., G. Bom, G. Rosborough, J. Marsh, and R. Williamson, Precise Orbit Computation 
for the Geosat Exact Repeat Mission, J. Geophys. Res., 95, C3, 2871-2885, 1990. 

Haines, B.J., G. Bom, J. Marsh, R. Williamson, A Summary of Precise Orbit Computation for 
the Geosat Exact Repeat Mission, John Hopkins APL, Technical Project, Vol. 10, No. 
4, 393-404, 1989. 

Hwang, C., Modeling Sea Surface Topography in Satellite Orbit Determination (abstract), EOS 
supplement, p. 89, April 23, 1991. 

Hwang, C., Modeling Sea Surface Topography in Satellite Orbit Determination (abstract), EOS, 
72, 17, 89. 


89 


Kim., J. and R.H. Rapp, The Development of the July 1989 1° x 1° and 30’ x 30' Terrestrial Mean 
Free- Air Anomaly Data Bases, Report No. 403, Dept, of Geodetic Science and 
Surveying, The Ohio State University, Columbus, January 1990. 

Koblinsky, C., J. Marsh, B. Beckley, A. Brenner, R. Williamson, Geosat Orbit Replacement 
Software for Altimeter Geophysical Data Records, GEM-T2 Ephemerides for 
November 1986 to November 1988, Geodynamics Branch, Goddard Space Flight 
Center, Greenbelt, MD, 2077 1 . 

Lerch, F., An Improvement Error Assessment for the GEM-T1 Gravitational Model, NASA Tech. 
Memo., TM 100713, 1988. 

Levitus, S., Climatological Atlas of the World Ocean, NOAA Professional Paper 13, U.S. Govt. 
Printing Office, Washington, D.C., 1982. 

Lisitzin, E., Sea Level Changes, Elsevier Scientific Publishing Co., New York, 1974. 

Marsh et al., Dynamic Sea Surface Topography, Gravity and Improved Orbit Accuracies From the 
Direct Evaluation of Seasat Altimeter Data, J. Geophys. Res., 95, C8, 13,129-13,150, 
August 15, 1990. 

Marsh J. et al.. The GEM-T2 Gravitational Model, J. Geophys. Res., 95,313, 22,043-22,072, 
1990. 

Mather, R.S., Determination of Some Dominant Parameters of the Global Dynamic Sea Surface 
Topography from GEOS-3 Altimetry, NASA Technical Memorandum 19558, Goddard 
Space Flight Center, Greenbelt, MD, 1978. 

Miller, L. and R. Cheney, Large Scale Meridionial Transport in the Tropical Pacific Ocean During 
the 1986-1987 El Nino from Geosat, J. Geophys. Res. 95, CIO, 17905-17,919, Oct. 
15, 1990. 

Nerem, S., B. Tapley, C.K. Shum, Determination of the Ocean Circulation Using Geosat 
Altimetry, J. Geophys Res., 95, C3, 3163-3179, 1990a. 

Nerem, R.S. et al., Determination of the Time- Variations of the Long Wavelength Dynamic 
Topography Using Geosat Altimetry, (abstract), EOS, 71, No. 43, p. 1266, Oct. 23, 
1990b. 

Pavlis, N.K., Modeling and Estimation of a Low Degree Geopotential Model from Terrestrial 
Gravity Data, Report No. 386, Dept, of Geodetic Science and Surveying, The Ohio 
State University, Columbus, p. 178, 1988. 

Pavlis, N.K. and R.H. Rapp, The Development of an Isostatic Gravitational Model to Degree 360 
and its Use in Global Gravity Modeling, Geophys. J. Int., 100, 369-378, 1990. 

Press, W.H., B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes, The Art of 
Scientific Computing. Cambridge University Press, 1986. 

Putney, B. et al., Earth Gravity Model Development at NASA/GSFC: Preliminary Results from 
GEM-T3S, EOS, p. 89, April 23, 1991. 

Rapp, R.H. and M. Kadir, A Preliminary Geoid for the State of Tennessee, Surveying and 
Mapping, 48, #4, 251-260, December 1988. 


90 



Rapp, R.H. and N.K. Pavlis, The Development and Analysis of Geopotential Coefficient Models 
to Spherical Harmonic Degree 360, J. Geophys. Res., 95, B13, 21,885-21,91 1, 1990. 

Rapp, R.H. et al., Consideration of Permanent Tidal Deformation in the Orbit Determination and 
Data Analysis for the Topex/Poseidon Mission, NASA Technical Memorandum 
100775, January 1991. 

Rapp, R.H., The Treatment of Permanent Tidal Effects in the Analysis of Satellite Altimeter Data 
for Sea Surface Topography, manuscripta geodaetica, 14, 368-372, 1989. 

Rapp, R.H., Y.M. Wang, N.K. Pavlis, Geoid Undulation Differences Between Geopotential 
Models, presentation at the meeting of the European Geophysical Society, Wiesbaden, 
Germany, April 1991a. 

Rapp, R.H., Y.M. Wang, N.K. Pavlis, Spatial and Spectral Differences Between Geopotential 
Models, EOS, (abstract), p. 89, April 1991b. 

Sacerdote, F. and F. Sanso, Spectral Calculus and Moving Average Operators on the Sphere.in 
Contributions to Geodetic Theory and Methodology, prepared by Section IV of the 
International Association of Geodesy, for the XX General Assembly of the IUGG, 
Vienna, 1991. 

Shibuya, K., Y. Fukuda, Y. Michida, Determination of Geoid Height at Breid Bay, East 
Antarctica, to appear in Journal of Geophysical Research, 1991. 

Torge, W., et al.. Long Range Geoid Control Through the European Traverse, Dtsch. Geod. 
Komm. Ser. B, 290, 1989. 

Wang, Y.M., Downward continuation of the free-air gravity anomaly to the ellipsoid using the 
grading solution Poisson’s integral and terrain connections, Rep. 393, Dept, of 
Geodetic Science and Surveying, The Ohio State University, Columbus, 1988. 

Wieser, M., The Global Digital Terrain Model TUG87, Internal Report on Set-up, Origin and 
Characteristics, Institute of Mathematical Geodesy, Technical University of Graz, 
Austria, 1987. 

Yi, Y. and R.H. Rapp, The October 1990 1° x 1° Mean Anomaly File Including an Analysis of 
Gravity Information from China, Internal Report, Dept, of Geodetic Science and 
Surveying, The Ohio State University, Columbus, February 1991. 


91 







Appendix A 

Start and Stop Times of the 76 Geosat Arcs Analyzed in This Study 


Arc 

# 

Start Time 
(UTC) 

Stop Time 
(UTC) 

No. of 
rr obs. 

rms 

(cm/sec) 

No. of * 
xovers 

rms* 

(cm) 

overlap 

rms(cm) 

1 

861 108 

861114 

7057 

.5054 

1999 

51.87 


2 

861113 

861119 

6214 

.5127 

2053 

68.48 

22.40 

3 

861118 

861124 

6334 

.5150 

2167 

49.40 

38.38 

4 

861123 

861124 

7236 

.5918 

2233 

57.98 

14.46 

5 

861128 

861204 

8216 

.5175 

2283 

48.15 

26.58 

6 

861202 

861207 2018 

7732 

.5318 

2146 

58.98 

13.49 

7 

861207 2100 

861214 

8037 

.4822 

2419 

53.65 

Bum 

8 

861213 

861219 

8203 

.4955 

2409 

50.37 

15.78 

9 

861218 

861224 

7798 

.5400 

2419 

51.13 

13.06 

10 

861223 

861229 

8570 

.5574 

2460 

50.30 

6,84 

11 

861228 

870103 

8138 

.4788 

3446 

48.57 

7.99 

12 

870102 

870107 2240 

9334 

.4955 

2594 

47.97 

8.43 

13 

870107 2323 

870114 

8663 

.5095 

2842 

47.27 

Bum 

14 

870113 

870119 

9642 

.4970 

2505 

47.34 

9.04 

15 

870118 

870124 

9182 

.5277 

2672 

44.28 

6.64 

16 

870123 

870129 

8677 

.4930 

2915 

49.45 

16.77 

17 

870128 

870203 

7496 

.5161 

2919 

5 6.65 .• 

22.97 

18 

870202 

870208 

7400 

.5012 

2567 

52.16 

31.66 

19 

870206 

870212 

8654 

.4975 

2744 

47.00 

8.13 

20 

870210 

870216 

8466 

.5087 

2729 

57.31 

17.30 

21 

870214 

870220 1021 

11329 

.5051 

2987 

55.49 

36.44 

22 

870220 1103 

8702226 

8516 

.5160 

2345 

43.72 

Bum 

23 

870225 

870303 

9302 

.4895 

2712 

50.28 

18.48 

24 

870302 

870308 

9163 

.5140 

2953 

55.99 

51.65 


68.40 
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25 

870307 

870313 

8674 

.5021 

2697 

44.71 

26 

870312 

870318 

8494 

.4963 

2784 

45.92 

27 

870317 

870323 

8817 

.5569 

2813 

51.26 

28 

870322 

870328 

7677 

.5679 

2898 

41.80 

29 

870327 

870402 

8920 

.5261 

2560 

43.21 

30 

870401 

870407 

9371 

.5336 

2348 

44.34 

31 

870406 

870412 

8675 

.5521 

2704 

43.77 

32 

870411 

870417 

9916 

.5169 

2687 

47.53 

33 

870416 

870422 

9393 

.5268 

2309 

44.51 

34 

870420 

870426 

10202 

.5378 

2106 

45.54 

35 

870424 

S70430 

11321 

.5371 

2224 

47.47 

36 

870428 

870503 1713 

9874 

.5502 

2165 

48.82 

37 

870503 1755 

870510 

12071 

.5036 

2732 

44.23 

38 

870509 

870515 

12247 

.5233 

2411 

47.59 

39 

870514 

870520 

12082 

.5133 

2314 

54.91 

40 

870519 

870525 

9322 

.5165 

2169 

46.04 

41 

870523 

870529 

9003 

.5005 

2041 

53.32 

42 

870527 

870602 

7394 

.4860 

2130 

55.02 

43 

870531 

870605 1829 

8359 

.5310 

2198 

56.52 

44 

870605 1912 

870612 

10995 

.5574 

2929 

47.67 

45 

870611 

870617 

10406 

.5367 

2112 

58.22 

46 

870616 

870622 1910 

11446 

.5626 

2840 

60.52 

47 

870622 1951 

870629 

9741 

.5172 

2386 

49.34 

48 

870628 

870704 

9663 

.5338 

2234 

51.30 

49 

870703 

870709 

10502 

.5268 

2045 

47.38 

50 

870708 

870714 

9588 

.4877 

1825 

51.12 

51 

870713 

870719 

9.166 

.5357 

1896 

52.17 


12.67 

19.00 

30.95 

23.28 

20.74 

24.20 

34.10 

32.38 

13.71 

9.39 

35.57 
Bum 
3.64 
48.88 

61.58 
3.34 
52.31 

62.07 
Bum 
15.02 
77.35 
Bum 
18.52 
27.94 

13.07 

50.07 
10.16 
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52 

870717 

870723 

10436 

.5115 

2009 

43.90 

11.20 

53 

870721 

870726 1230 

10547 

.5230 

1746 

43.76 

Bum 

54 

870726 1312 

870802 

11395 

.5245 

2697 

56.78 

59.59 

55 

870801 

870807 

10802 

.4795 

2478 

49.66 

33.98 

56 

870806 

870812 

10045 

.4920 

2491 

48.59 

64.94 

57 

870811 

870817 

11484 

.5146 

2525 

60.00 

77.38 

58 

870816 

870822 

9643 

.4886 

2383 

57.49 

49.40 

59 

870820 

870825 1517 

8930 

.5356 

2057 

50.25 

Bum 

60 

870825 1558 

870901 

10300 

.5258 

2274 

55.67 

81.86 

61 

870831 

870906 

9673 

.5661 

2259 

56.82 

74.11 

62 

870905 

870911 

9611 

.4992 

2220 

56.16 

81.92 

63 

870910 

870916 

10600 

.5189 

2203 

51.04 

60.55 

64 

870915 

870921 

9805 

.5456 

1877 

54.23 

82.27 

65 

870919 

870924 1614 

8752 

.5443 

1880 

64.12 

Bum 

66 

870924 1655 

871001 

10921 

.5257 

2479 

58.95 

59.84 

67 

870930 

871006 

10156 

.5407 

2142 

52.49 

48.60 

68 

871005 

871011 

9684 

.5402 

1890 

52.08 

44.78 

69 

871010 

871016 

11660 

* 

.5076 

1677 

62.86 

15.45 

70 

871015 

871021 

11870 

.5252 

1940 

60.49 

9.61 : 

71 

871019 

871024 1908 

10911 

.5340 

1810 

57.68 

Bum 

72 

871024 1951 

871031 

11628 

.5876 

2092 

57.25 ' 

8.56 ' 

73 

871030 

871105 

10153 

.5655 

2193 

65.18 

8.70 i 

74 

871104 

871110 

10137 

.5591 

2252 

58.24 

26.17 

75 

871108 

871114 

10855 

.5386 

1955 

58.64 

5.69 

76 

871112 

8711180017 

11131 

.5458 

1324 

63.90 



1 st year average 52.16 31.69 


Table from Koblinsky eL al (1990) 


94 





